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ABSTRACT 


Desirable  properties  of  linear  multistep  methods  (LMM)  can  be 
optimized  by  viewing  those  properties  as  functional  values  and  the  LMM 
possessing  those  properties  as  points  in  a  domain  space.  This  study 
conducts  two  such  optimizations  numerically.  The  first  is  a  search 
for  relatively  stable  explicit  LMM  and  the  second  is  a  search  for 
stiffly  stable  implicit  LMM.  Near-optimal ly  relatively  stable  explicit 
LMM  are  found  for  orders  four  through  nine. 

In  the  second  study,  the  concept  of  A(cc,  ^-stability  is 
introduced  for  stiffly  stable  LMM.  It  recognizes  the  need  for  large 
regions  of  absolute  stability  in  the  left  half  plane  and  the  need  for 
a  region  of  accuracy  about  the  origin  defined  by  the  region  of  relative 
stability.  An  economical  means  of  determining  the  region  of  relative 
stability  is  developed  and  used^_  Nearly-optimal  A(ct,  r)-stable 
implicit  LMM  are  found  for  orders  four  through  six  for  a  variety  of 
classes  determined  by  fixed  error  constants  C  +-j/a(l). 
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CHAPTER  I 


PRELIMINARIES 


1.1  Introduction 


Consider  the  initial  value  problem 


y'  =  f(x,  tj) ,  xc[a,  b],  ij{ a)  =  ijr 


(1.1-D 


where  if,  tj^,  and  f(x,  ij)  are  in  Rn  and  sufficient  conditions  are 
placed  on  f  to  ensure  a  unique  solution  exists.  We  seek  numerical 
methods  to  solve  the  problem  above  more  accurately  and  to  solve 
larger  classes  of  such  problems.  Especially  when  n  >  1  in  the 
problem  given,  it  can  be  roughly  classified  as  stiff  or  non-stiff. 

Vie  investigate  both  classes  and  find  improved  methods  for  solving 
members  of  each  class. 

Two  basic  approaches  to  a  numerical  approximation  of  the 
solution  to  (1.1-1)  are  the  linear  multistep  and  Runge-Kutta  methods. 
Most,  if  not  all,  popular  approaches  have  stemmed  from  one  or  both 
of  these  methods.  We  limit  the  ensuing  investigation  to  linear 
multistep  methods  (LMM). 


1.2  Linear  Multi  step  Approach 

When  it  is  not  possible  to  solve  the  continuous  system 
(1.1-1)  exactly  we  must  be  content  with  a  discrete  model  of  the 
given  system.  We  restrict  ourselves  to  an  evenly  spaced  mesh  of 
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the  interval  [a,  b].  Denote  this  mesh  by  {Xq  =  a,  •*•  , 


x-  =  a  +  ih,  •••  ,  xn  =  b}  so  that  we  have  constant  mesh-width 


b  -  a 


N  ‘ 


Define  </n  =  y(*n)  and  y ^  =  f{xn,  y  ),  and  denote  our  respective 


approximations  by  y^  and  y*  =  f(x  ,  y  ).  Our  discrete  model  then 

n  n  n  n 


consists  of  (y- }i  *  0,  ,  N}. 


A  linear  k-step  method  expresses  a  linear  relationship 


between  yn+-|>  y^+-j  and  the  previous  k  values  of  y^  and  y!.  This 
linear  relationship  takes  the  form 


k-1  k-1 

Vi  =  T0  Vn-i +  h  T„,  W-i 


(1.2-D 


for  real  values  of  a.  and  ,  where  |  +  i  S^  -j  |  f  0.  Returning 


to  the  continuous  model,  consider 


i\—  i  rv  • 

l  a.ij  .  +  l  &■!/'  . 
i=_l  iJn-i  >-j  ^Vn-i 


where  a_-j  =  1.  If  «/(x)  is  sufficiently  differentiable,  by  a  Taylor 

expansion  about  x  for  ij  .  and  u‘  . ,  i  =  -1,  0,  1,  •••  ,  k-1 ,  v/e 
n  Jn-i  ^n-i 


J./Vn-i  +  -  l0  (1.2-2) 


where  the  constants  C.  are  determined  from  the  coefficients  a.  and  3^ 
of  the  LMM  (1.2-1).  We  have  the  following  expressions  for  the  : 


k-1 


c° "  ill v 


k-1 


k-1 


Cl  =  I  (-l)a,  +  I  B- : 
1  i=-l  1  i=-l  1 


(1.2-3) 


CP  =  h  tL,  H)P°i  +  Tp^  j.,  for  P  x  2>  3’  "• 


The  order  of  a  LMM  (1.2-1)  is  defined  to  be  the  smallest 
integer  p  such  that  Cj  =  0  for  0  _<  i  _<  p  and  C  +1  t  0.  When  the 
order  is  at  least  one  the  LMM  is  said  to  be  consistent.  The  following 
concept  of  convergence  is  critical  for  any  useful  LMM. 

Defini ti on:  A  LMM  is  said  to  be  convergent  if  for  any  problem 
(1.1-1)  for  which  we  are  guaranteed  a  unique  solution,  we  have 


lim  y  =  y(x) 
n-*0  n 
nh  =  x-a 


for  all  xe[a,  b]  and  for  all  solutions  {yn>  of  (1.2-1)  satisfying 


starting  conditions  y.  =  n^(h)  where 


lim  n.(h)  =  yn,  i  =  0,  1,  ,  k-1 

h-K>  1  u 


With  the  LMM  (1.2-1)  we  associate  two  polynomials  p(z)  and 
o(z).  Define 
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p(z)  =  l  a.z 
i=-l  1 


k-(i+l) 


(1.2-4) 


k-(i+l) 


o(z)  =  l  0- z 
i=-l  1 


where  a  ^  =  -1.  Then  any  L 1414  (1.2-1)  uniquely  determines  both  p  and 
o  and  conversely.  If  p(z)  has  no  roots  greater  than  one  in  modulus 
and  if  all  roots  of  modulus  one  are  simple,  the  associated  LMM  is 
said  to  be  zero-stable.  We  have  the  following  theorem  (stated 
without  proof)  from  Dahlquist  [5]. 

Theorem  1.2-1.  A  LMM  is  convergent  if  and  only  if  it  is 
consistent  and  zero-stable. 

Formula  (1.2-1)  provides  an  explicit  naans  of  obtaining 
yn+1  when  0  -j  =  0.  The  resulting  LMM  are  called  predictors.  When 
0  ^  f  0  the  formula  defines  yR+^  implicitly.  In  this  case  iterative 


val ues 


(yMr 


are  computed  via  the  relationship 


v[s+l]  _  w  [s], 
yn+l  “  ^Vl' 


(1.2-5) 


where  <j>  is  defined  by  the  right-hand  side  of  (1.2-1).  As  long  as  <•> 


is  contractive  sequence 
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(y 


S=1 


will  converge  to  y  This  iteration,  normally  used  to  obtain  the 
value  of  y^  when  £_-j  j*  0,  prompts  the  name  correctors  for  this 
class  of  LMM.  For  any  k-step  method  k  starting  values  are  required 
to  implement  the  recursive  calculation  of  the  y.. ,  i  _>  k.  These 
values  are  normally  found  by  some  explicit  LMM  or  R-k  method. 

Denote  the  class  of  n**1  order,  m  step  predictors  and 

correctors  respectively  by  P{n,  m)  and  C(n,  m) .  For  example  the 
th 

n141  order  Adams-Moulton  formula  is  a  member  of  c{n,  n-1).  Similarly 
tfie  n^1  order  Adams-Bashforth  formula  is  a  mender  of  p(n,  n).  We 
search  within  p(n,  n),  n  =  4,  •••  ,  9  for  imp roved  methods  on 
non-sti  ff  problems  and  within  C(n,  n) ,  n  =  4,  5,  6  for  improved  methods 
on  stiff  problems. 

1 . 3  Coefficient  Matrix  Derivation 

Let  l  be  a  LMM  of  the  class  C(n,  n).  As  in  (1.2-1)  we 
normally  assume  a =  1  to  eliminate  ambiguity  between  equivalent 
LMM.  Then  from  (1.2-3)  we  have  n+1  independent  equations  relating 
2n  +  1  parameters.  This  leaves  n  unspecifie-d  (free)  parameters  to 
determine  L.  We  take  these  free  parameters  to  be  a-j ,  •••  ,  an  -j, 
and  £_^.  In  (1.2-3)  Cq  =  0  determines  Oq  and  C-j  =  ••*  =  Cn  =  0 
yields  n  x  n  matrices  B  and  D  and  an  n  x  (n+1)  matrix  A  such  that 
DA[1,  a-j,  •••  ,  an_i>  =  B[Bq>  •••  ,  Bn_-|]T.  These  matrices 

have  the  form 


1 


A  = 


-1 


(-l)n_1  1  2n 


(n-1) 

(n-1)2 


(n-1 )n 


-1 

2 


(-D 


n 

n 


(1.3-1) 


1  1  1 
0  1  2 


B  = 


0 


1 


1 

(n-1) 

« 

(n-1)"-’ 


and  D  =  (djj)  is  diagonal  with  d.^  =  1/i.  In  this  section  we  define 
B_1  recursively,  that  is,  we  define  the  B-1  for  C(n,  n)  from  informa¬ 
tion  concerning  the  B"1  for  C(n-1,  n-1).  Thereby  we  have  the  means 
to  calcuate  the  exact  rational  entries  of  the  n  x  (n-:l)  matrix 
B~^DA  and  so  all  coefficients  of  L  are  explicitly  determined  by  the 
n  free  parameters. 

Define  the  n  x  n  matrix  Sn  recursively  as  follows.  Let 
S1  =  (1),  S2  =  [o  l), 

and  Sn  =  (c1?-)  where 
•  3 

°1  1  =  (n_^!  »  °i  -j  =  0  for  2  _<  i  _<  n. 


o"  „  =  1  for  1  < 
i,  n  - 


<  n. 


j 


j-1 


for  2  <  j  <  n-1 , 


l 

i 


: 


t 

t 


i 


\ 

i 

§ 


i  m 


nwt*  ,****■*< 
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and  o'?  .  =  a"'1  -  ,  +  (n-1)  o1?"^ .  for  1  <  i  <  n-1  and  2  <  j  <  n-1. 
i,  J  i,  J-l  i ,  J  -  -  —  — 

Denote  the  for  C{n,  n)  by  B-^  =  (b1?.).  Then  we  have 

n  1 3 


„  (-Df+j  o1:-  . 
bi,  j  =  TTTTTtlPTjr1  f»r2<1.  j<n, 


bi,  l  ‘  ’•  b?,  o 


n 


i 


n 


bV  .  for  2  <  j  <  n, 
Z=2  *»  3  ~  ~ 


and 


b.  1=0  for  2  <  i  <  n. 
1>1  —  — 


-1. 


When  computing  the  exact  rational  entries  of  the  matrix  DA  for 
higher  orders  it  is  necessary  to  increase  the  machine  precision  and 
use  careful  programing.  He  have  calculated  rational  entries  for 
these  matrices  of  coefficients  through  order  18.  In  Appendix  A  we 
list  those  matrices  for  orders  two  through  nine. 


CHAPTER  II 


STABILITY  ANALYSIS 
2.1  A  General  Consideration 

Many  questions  concerning  the  stability  of  LKM  reduce  tc 
determining  whether  certain  polynomials  have  roots  less  than  one 
in  modulus.  Much  effort  towards  this  end  has  been  devoted  to 
transformed  polynomials  whose  roots  lie  in  the  open  left  half  plane. 

Ke  say  a  polynomial  is  of  the  Schur  type  if  it  has  roots  less  th3r. 
one  in  modulus  and  of  the  Hurwitz  type  if  it  has  all  roots  lying  in 
the  open  left  half  plane.  Several  forms  of  necessory  3nd  sufficient 
conditions  exist  in  the  literature  for  a  polynomial  tc-  be  Hurwitzian 
[23,  26].  Translation  to  necessary  and  s  ufficient  conditions  for 
Schur  polynomials  normally  leads  to  intractable  criteria.  Ke  nanaoe 
to  find  one  exception  to  this,  although  in  the  general  case  it 
provides  sufficiency  only.  First  we  establish  several  preliminary 
results. 

fi 

Lensna  2.1-1.  Let  A  be  an  (n+l)x  n  matrix  and  CcR" .  Also 

let  the  set  of  points  xcP.n  satisfying  the  system  or  linear  inequalities 

Ax  >  C  define  an  n-dimensional  simplex  £  with  non-emoty  interior. 

—  n 

Let  t  be  the  set  of  all  points  i/eR0  such  that  ky  <  C.  Then  t  is 
ercpty. 

Proof:  Let  Xg  be  a  point  interior  to  En  and  suppose  there 
exists  a  point  £fgCT.  Consider  the  line  segment  aXg  +  (l-a)^  where 


9 

0  <o  <  1.  This  line  segment  must  intersect  each  of  the  n+1 
hyperplanes  defining  En  since  AxQ  >  C  and  A tjQ  <  C. 

Now  extend  this  line  segment  to  a  >  1.  Since  En  is  bounded 
this  extension  must  intersect  at  least  one  of  the  defining  hyperplanes 
a  second  time  at  a  point  distinct  from  the  first  point  of  intersec¬ 
tion.  Thus  we  have  the  line  axQ  +  (l-a)f/0  contained  in  one  of  the 
defining  hyperplanes  of  Zn«  This  contradicts  the  assumptions  that 
AxQ  >  C  and  A tjQ  <  C.  QED. 

Lemma  2.1-2.  Let 

P(z)  =  I  B.z1 
i=0  1 

be  a  Hurwitz  polynomial.  Then  all  >  i  =  0,  1,  ,  n  have  the  same 

sign. 

Proof:  If  all  roots  of  P  have  negative  real  part,  then 

iDdsl 

P(z)  =  B  H  (z  +  r.)  n  ((z  +  C,)2  +  d  2) 
n  i=l  1  i=l  11 

where  r^  >  0  and  >0  so  that  all  coefficients  take  the  sign 
of  3n*  QED. 

Lemma  2.1-3.  Consider  the  mapping 


which  maps  the  unit  disk  onto  the  left  half  plane.  For  any  polynomial 
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’  -*r ' 


define  a  new  polynomial 


P(z)  =  (z-l)n  p[gf. 


Then  for  the  polynomial  transformation  T  :  p  -*•  P  we  have  T  = 
where  I  is  the  identity  transformation. 


Proof: 


?(p)  =  T(P)  =  (z-l)n  p[|±|]  =  (2-l)n(^§T]n  pjf2- 


=  2n  p(z)  =  2n  I(p).  QED. 


Consider  the  real  monic  polynomial 


i 

p(z)  =  l  a.z  ,  a =  1. 
i=0  1  n 


& 


Let  the  coefficient  vector  of  p  be  a  =  (cxq,  >  an  -|)  and  «j  be 

the  coefficient  vector  of  (z-l)n~J  ( z+1 ) ^ .  Define  Sn  to  be  the 

n-dimensional  simplex  determined  by  a.,  j  =  0,  1,  •••  ,  n. 

J 

Theorem  2.1-4.  The  coefficient  vector  of  any  nth  order  real 
monic  Schur  polynomial  lies  in  S  . 

Proof:  Let 


i 

p(z)  =  l  a.z  ,  a  =  1, 
i=0  1  n 


be  a  real  monic  Schur  polynomial.  Transform  the  roots  of  p  to 
the  left  half  plane  via  the  mapping 


OS'- 


l  I 


to 


Define  the  new  polynomial 


2  +  1 
z  -  1  ' 


p(z)  =  I  e,?.1 
i=0  1 


as  in  Lemma  2.1-2.  Then 


P(z)  -  I  a,.(z-l)"'1(z+l)' 
i=0  1 


n-j  fn-j'l  f  j 


6-  =  l  a.  I  kj  [i-kj(-l)n"j“k,  i  =  0,  1,  ,  n. 

1  j=0  J  k=0 

0<i-k<j 


From  this  expression  denote  the  coefficient  of  z1  in  (z-1  )n-J(z+l )' 

by  and  let  Mn  =  (m^),  0  _<  i,  j  £  n.  We  have  then  a  matrix 

T  T 

equation  Mn(a0,  ,  an)  =  (Bq,  •••  >  $n)  .  We  note  that  Mn  can 
also  be  described  by 


=  bJ>  i  =  °»  1 »  **•  >  n 


"nr1-  j  =  0>1 


n  n  ,  n  ,  n  4  —  n  i  i 

mij  -  "i,  o-1  mi+l  >  j-1  +  Vl.  J  •  ’  -  “•  ’•  "•  ’  n'1 

j=l,2,  •••  ,n. 


.•j(cC  V>vV»jy,'  Jv  k  ■£>  >&  v.>  v 


Lensiia  2.1-3  provisos  M1"  ~  2ni  . ,  5  w hare  I  . ,  is  the  (nil)  x  (m-1) 
r  n  ni  1  n-i-1  ' 

identity  matrix. 

Since  p-'z)  is  a  real  manic  polynomial ,  the  first  n  elements 

J.  L. 

in  the  j  n  column  of  M  'is  the  vector  a.  described  in  the  discussion 

n  j 

previous  to  the  statement  of  the  tiieoiem.  Thus  taking  the  product 
2  p 

M  corresponds  to  evaluative]  rvl  hype rpl sues  in  R'  at  the  corff  cient 

2  r. 

vectors  «•>  j  -  0,  1,  •  ,  n.  Since  H‘  =  2‘In+[  each  of  the 

coefficient  vectors  a.  lies  on  all  but  one  of  the  r.+l  hype rpl arcs . 

J 

They  therefore  describe  S  . 

J  n 

Let  f.  represent  the  k  hyperplane  where  for  xcRn  the  k1'*1 

l\ 

T 

row  of  M  (x, ,  5  x  ,  1)'  defines  f.(x).  We  define  the  positive 

iii  n  k 

side  of  the  k  "  hyperplane  to  consist  of  those  points  xcR  for  which 

2  q 

f.  (x)  >  0.  =  2' J  ,,  iinolies  each  a.  lies  on  the  positive  side 

k  n  n-H  ■  j 

of  its  opposite  face.  Thus  for  each  xeRn  wo  iiave  xcSn  if  and  only 
if  f.  fx)  is  positive  for  each  k  =  1,  >  ri  +  1. 

Notice  if  0  is  the  zero  u*  Rn  v/e  have  ^(6)  >  0  for  each 
k  -  1,  ,  n  +  1 ,  so  that  continuity  of  the  fj,  provide  Sn  with  a 

non-empty  interior. 

Since  p  is  a  Schur  pclynonn'al,  P(z)  is  a  llurwitz  polynomial 
and  lemrna  2.1-2  provides  that  all  g^  are  of  the  same  sign.  Lemma 
2.1-1  excludes  the  possibility  that  all  g^  are  negative.  Since 
B.j  =  (aQ,  ,  •••  ,  en_])  we  ^ave  the  coefficient  vector  of  p(z) 
in  Sn<  QED. 

This  Theorem  is  in  fact  necessary  and  sufficient  when  n  _<  2. 
We  obtain  the  following  useful  corollary. 


Lei. 


Corollary  2.1-5. 


pU) 


=  V 
1=0 


and  M  be  as  in  the  proof  of  Thoc •  e:::  0.1-4.  Let 
n 


P( 


•'2) 


n 

X  cy 

i=o  1 


where  M  (ar.  a,  >  ,  a  )  =  (P0»  (S15  *’•  5  3  )‘.  Then  if  p(z) 

is  a  Sehnr  polynomial  ail  the  have  the  same  sign 

Proof:  From  the  proof  of  the  Theorem  wo  see  that  P(z)  is  a 
Hunvitz  polynomial.  The  result,  follows  upon  application  of  Lemma 
2.1-2.  QEI). 


2 * 2  Stability  fc r  Nun-Sti ff  LMM 

In  Chapter  I  wo  mentioned  convenience  of  a  LMM  and  its 
equivalence  to  consistency  and  zero-stability.  There  we  were 
concerned  about  what  happens  to  the  error  through  successive 
calculations  at  a  fixed  point  within  the  interval  of  integration  as 
h  ->  0  and  n  •*  ».  The  definition  of  convergence  requires  this  error 
to  go  to  zero  on  all  problems  of  a  certain  class  giver,  sufficiently 
accurate  starting  values.  Now  we  concern  ourselves  with  what  happens 
to  the  error  through  successive  calculations  using  a  fixed  step 
length  h  as  we  proceed  through  the  interval  of  integration.  Odeh 
and  Liniger  [30]  refer  to  this  concept  as  fixed-h  stability. 

Defini ti on:  A  LMM  is  called  fixed-h  stable  if  the 
accumulated  truncation  error  in  solving  the  model  equation  ij'  =  Xij , 
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where  X  is  a  complex  constant,  remains  bounded  as  n  ->  °°  for  a 
non-vanishing  range  of  values  of  the  constant  q  ~  Ah. 

Here  we  come  to  a  disjunction  in  the  stability  analyses  for 
stiff  and  non-stiff  LMM.  Non-stiff  methods  normally  have  bounded 
regions  of  fixed-h  stability  in  the  q-plane  and  stiff  methods  require 
those  regions  to  be  unbounded.  We  pursue  the  analysis  for  stiff 
methods  more  fully  in  the  next  section. 

The  concepts  of  absolute  stability  and  relative  stability 
are  more  commonly  used.  We  denote  the  characteristic  polynomial 
of  a  LMM  by  ll(z,  q)  =  p(z)  +  qa(z)  where  p  and  a  are  taken  from 
(1.2-4)  and  q  is  defined  as  in  the  previous  definition  from  the 
model  test  equation  ij'  =  Xij. 

Defi ni ti on :  A  LMM  is  absolutely  stable  at  a  given  point  q^ 
if  all  roots  of  n(z,  q^)  lie  inside  the  unit  circle.  The  region  of 
absolute  stability  for  a  LMM  is  the  set  of  all  points  qQ  in  the 
q-plane  at  which  the  method  is  absolutely  stable. 

Absolute  stability  forces  solution  of  the  model  test 
equation  to  have  a  decreasing  global  error  and  is  therefore  more 
restrictive  than  fixed-h  stability. 

For  any  consistent  and  zero-stable  method  L ,  n(z,  0)  =  p(z) 
has  a  simple  root  at  z  =  1.  This  root  is  called  the  principal 
root  and  is  denoted  by  r^.  As  long  as  the  leading  coefficient  of  a 
polynomial  is  non-zero,  its  roots  are  continuous  functions  of  the 
coefficients.  Thus  we  may  follow  this  principal  root  r^  to  non-zero 
values  of  q.  For  q  sufficiently  near  zero  we  have 
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q  ,  m  K+l, 

r1  =  eM  +  0(q  ) 


(2.1-1) 


whenever  l  is  of  order  k.  For  an  argument  establishing  this 


relationship  see  Lambert  [27,  p.  66], 


Definition:  A  LMM  is  relatively  stable  at  a  given  point 


if  the  root  of  n(z,  q^)  of  largest  magnitude  is  the  principal  root. 


The  region  of  relative  stability  for  a  LMM  is  the  set  of  all  points 


q^  in  the  q- plane  for  v/hich  the  method  is  relatively  stable. 


There  are  numerous  definitions  of  relative  stability  given 


in  the  literature  (for  example  see  [27,  p.  68]);  in  light  of  (2.1-1) 


the  definition  above  is  very  useful.  We  take  a  closer  look  at 


alternate  definitions  in  Chapter  IV.  Relative  stability  limits 


the  rate  of  growth  of  the  global  error  to  approximately  that  of  the 


solution.  We  will  examine  this  relative  error  more  closely  in 


Chapter  III. 


For  methods  which  are  consistent  and  zero-stable  there  will 


not  be  regions  of  absolute  stability  near  the  origin  in  the  right 


half  plane  as  there  are  with  relative  stability,  while  in  the 


left  half  plane  the  regions  of  absolute  stability  are  larger  than 


those  of  relative  stability.  These  conclusions  follow  since  for  q 


arbitrarily  near  the  origin  in  the  right  half  plane  we  have  r^  >  1 , 
and  for  q  in  the  left  half  plane  r-j  <  1.  When  solving  a  problem  with 
q  lying  in  the  left  half  plane  and  inside  the  region  of  absolute 


stability  but  outside  the  region  of  relative  stability,  the  global 


error  goes  to  zero  but  possibly  much  slower  than  the  solution 


itself.  We  illustrate  this  possibility  with  an  example  in  Chapter  IV. 
Our  results  must  then  be  viewed  with  this  is  mind. 


Throughout  the  remainder  of  the  section  we  assume  p(z)  and 
o(z)  have  no  roots  in  common  and  that  a(z)  has  no  roots  of  unit 
magnitude.  To  find  the  boundary  of  the  region  of  absolute  stability 
we  take  the  image  of  the  unit  circle  under  the  map 


z  -*■  q 


This  boundary  is  not  always  a  simple  closed  curve.  The  boundary  of 
the  region  of  relative  stability  is  found  numerically  by  tracking 
the  principal  root  on  various  rays  emanating  from  the  origin.  These 
notions  are  considered  more  carefully  in  Chapter  IV.  The  following 
result  is  often  used,  many  times  casually,  and  we  now  state  it 
formally. 

Lemma  2.2-1.  The  boundary  of  the  regions  of  absolute  and 
relative  stability  are  symmetric  about  the  real  axis. 

Proof:  If  q  and  q  are  complex  conjugates  we  have  H(z,"  q)  = 
n(z,  q)  so  that  zQ  is  a  root  of  Jl(z,  q)  if  and  only  if  zQ  is  a  root 
of  n(z,  q).  nus  the  sets  of  moduli  of  the  roots  of  H(z,  q)  and 
n(z,  q)  correspond  and  also  |z^(q)|  =  |z^(q)|  where  z-j(q)  is  the 
principal  root  of  n(z,  q).  Since  the  stability  boundaries  are 
based  on  the  moduli  of  the  roots  of  n(z,  q)  the  result  follows. 

QED. 

We  now  take  a  closer  look  at  the  boundary  of  the  region  of 
absolute  stability  and  develop  several  results.  We  use  the  following 
not.ati  on. 


p(eieV  Pi(0)°r(Q)  "  Pr(6)oi(9) 
ifl  '  =  ?  5 

k  o(e  8)J  °p  (0)  +  o/ (0) 


we  can  express  a.  as 

D 


_1  fp  (0)a  (0)  +  p.(0)a-(0)' 
ctA  =  cot  I  .  /  n\~  T7T\  ~~rzvz  rTTT 


pjJeTaT[e)  -  atJeOpi(oTj‘ 


(2.1-2) 


Lemma  2,2-2.  For  any  convergent  LMM  there  is  a  deleted 
neighborhood  of  zero  in  which  ^ (9)  f  0. 

Proof:  We  have 


p!(e)  =  l  (k-j-l)a.  cos  (k-j-l)e. 
1  j=-l  J 


pUO)  =  l  (k-M)ot.  =  p'(l)  f 
1  j=-l  J 


since  the  LMM  is  convergent  and  therefore  has  a  simple  root  at 
z  =  1.  Since  p.!  is  continuous,  there  is  a  neighborhood  Ng  of  zero 
in  which  pi  f  0.  Thus  for  SeNg,  9^0,  we  have  p^(0)  ?  0.  QED. 
Lemma  2.2-3.  For  any  consistent  LMM 


pr(e) 

lim  - — tvt  -  0. 
9-K)  pi 


mum 
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Proof: 


Pr(e) 

1  im  — rz~v 

e->0  V9' 


i=_i  o 
=  lim  4 


k-1 

l  a.  cos  (k-j-1 )e 


0->0 


k-1 

I  a.  sin  (k-j-1 )9 
j=-l  J 


k-1 


l  (j+l-k)a,  sin  (k-j-1 )0 

iz-j _ 3 _ 

=  lim  k-1 

e-’-O  (k-j-1  )a.  cos  (k-j-1  )G 

j=-l  3 


since  consistency  insures 


k-1 

l  (k-j-1 )a.  /  0. 

,i*-l  3 


QED. 


Theorem  2.2-4.  For  any  convergent  LMM 


lim  a  =  n/2. 
6+0  9 


Proof: 


lim  a 

e->o 


e 


lim  cot 
0+0 


p  a  +  *p.o. 
r  r  ii 


pa.-  p.o 
r  i  i  r 


lim 

6+0 


(p  /p.)a  +  a. 
*  r  Ki  '  r  i 

<pr/ci)oi  -  °r 


=  n/2. 


since  o(l)  f  0.  QED. 
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Thus  we  see  che  boundary  locus  of  the  region  of  absolute 
stability  always  leaves  the  origin  along  the  imaginary  axis.  Indeed 
then,  if  the  non-principal  roots  of  p  are  strictly  inside  the  unit 
circle,  there  will  be  a  non-vanishing  interval  of  absolute  stability 
to  the  left  of  the  origin  and  there  will  be  a  non-null  region  of 


absolute  stability  in  the  left  half  plane. 

For  0  =  0,  n,  since  p  and  o  are  real  polynomials,  the  boundary 
locus  points  will  be  real.  At  these  values  of  6  we  can  compute  the 
boundary  locus  points  in  general. 

Theorem  2.2-5.  The  boundary  locus  points  at  6  =  0,  H  for  the 
region  of  absolute  stability  of  any  convergent  LMM  are  the  origin 
and  p(-l)/c(-l),  respectively. 

Proof:  For  e  =  0,  convergence  provides  a  simple  root  of 
p(z)  at  z  =  1.  Consistency  implies  p'(l)  =  o(l),  thus  o(l)  =  0 
forces  a  multiple  root  of  p  at  z  =  1.  We  conclude  the  origin  is  the 
boundary  locus  point  at  8  =  0. 

For  0  =  II  we  have 


lim  Re 
e-n 


,  i0,l  p  a  +  p.o.  p  (n)a  (Jl)  ,  , , 

p(e  )  _  ,.m  r  r  i  _  Krv  rl  p(-1) 
- rrr~  ~  nm - * - n—  -  - o - — rr* 


a(e  )j  e+n  op  +  oi 


[ar(n)r 


apT) 


There  may  be  other  values  of  0  for  which  the  boundary  locus 
points  are  real.  For  those  LMM  of  interest  to  the  investigation 
reported  in  Chapters  IV  and  V  v/e  characterize  real  boundary  points 
in  the  theorem  following.  First  we  define  a  desired  property. 
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Definition:  A  LMM  is  strictly  zero-stable  if  the  associated 
polynomial  p(z)  has  a  -inple  root  at  z  =  1  and  all  remaining  roots 
interior  to  the  unit  circle. 

This  definition  is  clearly  stronger  than  the  zero-stability 
defined  in  Chapter  I. 

Theorem  2.2-6.  Let  L  be  a  strictly  zero-stable  LMM  for  which 

no  roots  of  a(z)  have  unit  magnitude.  Then  the  boundary'  locus  points 

of  the  region  of  absolute  stability  are  real  if  and  only  if  p(e*9)  and 
i  9 

o(e  )  lie  on  the  same  line  through  the  origin. 

Proof:  When  0=0  the  result  is  clear  since  both  p(l)  and 


*  : 


o(l)  are  real. 

In  the  remainder  of  the  proof  we  assume  8/0.  The  boundary 
locus  points  are  real  t 


( 


Im 


P(ei0) 


l  o(ei6)J 


=  0 


p .a  -  p  o. 
r  ti 

2  T  2 

°r  +  0i 


0-  P^r=  P^- 


Suppose  p. (6)  =0.  If  the  boundary*  locus  points  are  real 

then  Pr(e)o.(B)  =  0.  Since  6/0  and  L  is  strictly  zero-stable 

Pr(6)  /  0  so  that  0^(0)  =  0.  Then  both  p(el9)  and  o(e^9)  are  real  and 

i  6 

lie  on  the  same  line  through  the  origin.  Conversely,  if  both  p(e  ) 
i  9 

and  o(e  )  are  on  the  same  line  through  the  origin  and  p^e)  =  0, 

18  18 

both  p(e  )  and  o(e  )  are  real  and  the  boundary  locus  point  is  real. 

Suppose  0^(9)  =0.  If  the  boundary  locus  points  are  real 
then  p.(9)or(9)  =  0.  Since  or(0)  =  0  would  yield  a  root  of  o 
with  unit  magnitude  we  conclude  p.(e)  =  0.  Thus  both  p(el9}  and 


'll 
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a{elG)  are  real  and  lie  on  the  same  line  through  the  origin. 

Conversely,  if  both  p{e  )  and  o(e  )  lie  on  the  same  line  through 
the  origin  and  a  -  ( 8)  =  0  then  p.(9)  =  0  and  both  p(e16)  and  o(e19)  are 
real  resulting  in  a  real  boundary  locus  point. 

If  both  p-(G)  and  { 6)  are  nor.- zero  by  the  remarks  at  the 
beginning  of  the  proof  we  have  the  boundary  locus  point  at  6  real  if 
and  only  if 

Pr(0)  cr(9) 

pTTeT  =  oTTel  ‘ 

These  are  equivalent  to  p(elS)  and  o(el0)  lying  cn  the  same  line 
through  the  origin.  QED. 

2.3  Stability  for  Stiff  L,'»M 

When  n  >  1  in  problem  (1.1-1)  the  eigenvalues  of  the  system  can 
be  of  vastly  different  magnitudes  within  the  interval  of  integration. 
When  this  situation  occurs  stability  constraints  on  the  maximum 
stepsize  permitted  may  be  dictated  throughout  the  interval  by 
eigenvalues  whose  contribution  to  the  solution  of  the  system  becomes 
negligible  after  a  time.  Often  tie  stability  constraint  limits  the 
stepsize  so  severely  that  the  roundoff  errors  and  amputation  time 
involved  are  overwhelming.  It  would  be  nice  v/hen  this  occurs  to  have 
a  region  of  absolute  stability  extending  to  infinity  to  allow  us  to 
determine  stepsize  on  the  basis  of  accuracy  requirements  alone. 

This  is  the  essence  of  the  stiff  problem  and  its  solution. 


We  wish  then  to  find  LMMwith  infinite  regions  of  absolute 
stability  in  the  left  half  plane  without  sacrificing  accuracy  around 
the  origin.  First  we  review  the  more  important  definitions  from  the 
literature  characterizing  properties  of  such  methods.  One  such 
definition  has  been  provided  by  Gear  [9]  and  descrit-s  stiffly  stable 
LMM. 

Defi ni ti on:  A  LMM  is  stiffly  stable  if  in  the  region  R-j 
( Re q  <:  D)  it  is  absolutely  stable  and  in  R£  (D<Re  q  <  a,  |Imqj  <  G) 
it  is  accurate. 

The  region  of  accuracy  referred  to  here  is  not  precisely 
defined  and  therefore  this  definition  is  not  entirely  workable. 

Other  definitions  focus  only  on  the  infinite  region  of  absolute 
stability  and  omit  any  measure  for  a  region  of  accuracy  about  the 
origin.  The  following  definition  provides  a  generally  unobtainable 
standard  to  which  methods  may  be  compared.  It  was  introduced  by 
Dahlquist  in  1963. 

Defini tion:  A  LMM  is  called  A-stable  if  the  region  of 
absolute  stability  includes  the  open  left  half  plane. 

We  follow  here  with  a  somewhat  discouraging  result  also  due 
to  Dahlquist  [6]. 

Theorem  2.3-1.  The  order  of  an  A-stable  LMM  cannot  exceed 

ty/o. 

This  result  naturally  led  to  less  severe  definitions. 

Widlund  [41]  introduced  the  concept  of  A(a)-stability. 


Definition:  A  LMM  is  A(a)-stable,  ae(0,  11/2)  if  the  region  of 
absolute  stability  includes  the  wedge 


=  {z|  |Arg(-z) |  <  a,  z  t  0). 


(Here  Arg  zc.(-Jl,  n]V  zeC). 

A  method  is  A(n/2)-stable  if  it  is  A(a)-stable  for  all 
ae(0,  II/2)  and  A(0)~stable  if  it  is  A(a)-stable  for  some  (sufficiently 
small)  ac(0,  FI/2) . 

Cryer  [4]  completed  the  spectrum  of  A(ct) -stability  with  his 
definition  of  Ag-stabil i ty. 

Definition:  A  LMM  is  A^-stable  if  its  region  of  absolute 
stability  includes  the  negative  real  axis. 

Clearly  then  for  systems  with  real  eigenvalues  we  need 
Ag-stability  and  for  complex  eigenvalues  we  need  A(«)-stability  where 
the  a  is  determined  by  the  eigenvalues  of  the  system  we  are  solving. 
Cryer  [4]  gives  an  example  of  a  LMM  which  is  AQ-stable  but  not  A(0)- 
stable.  There  exists  then  a  series  of  strict  inclusions  for  the 
stability  classes  of  LMM  from  A-stability  through  Ag-stability. 

Cryer  [4]  shows  that  the  Adams- 1*10111  ton  formulae  of  order 
k  _>  2  are  not  A^-stable  but  he  does  show  that  there  exists  A^-stable 
methods  of  arbitrarily  high  order.  The  backward  differentiation 
formulae  implemented  by  Gear  are  not  stiffly  stable  for  order  k  _>  7. 
Our  investigation  of  LMM  for  use  on  stiff  systems  is  restricted  to 
implicit  methods  as  we  shall  see  presently.  First  a  variation  of 
a  result  from  Rodabaugh  [34]. 

Lemma  2.3.2.  Let  (a),  i  =  -1,  m  be  continuous 
functions  of  a  on  the  interval  (-«>,  b)  with 


6, 


-.-M- 


i  m 
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1  i in  X  ,(a)  =  0. 
a ->•-«> 

Assume  that  for  some  i, 

Tim  A..(a)  /  0. 
a->-°° 

Then  for  any  x  >  0  there  exists  an  N  such  that,  if  a  <  N  then 

m  m  . 

I  A.fajr'"-1  =  0  (2.3-1) 

i=-l  1 

has  a  root  with  absolute  value  larger  than  x. 

Proof:  Let  j  =  min  {i  |  i  _>  0  and  1  i in  A. (a)  f  0).  The  roots 

a->-«>  1 

of  (2.3-1)  are  identical  with  those  of 

rnH'1  +  l  ( X. ( a)/X  , (a) ) rm-1  =  0  (2.3-2) 

i=0  1 

provided  A_ -j ( a )  f  0.  Since  (2.3-2)  is  monic,  its  coefficients 

(A^(a)/A_^(a))  can  be  expressed  by  the  symmetric  polynomials  as 

polynomial  functions  of  its  roots.  However,  as  a  tends  to  negative 

infinity,  A.(a)/A  .(a)  increases  without  bound  so  at  least  one  of  the 
3  * 

roots  is  unbounded.  QED. 

Theorem  2,3-3.  Any  explicit  LMM  is  not  A^-stable. 

Proof:  Let  Ii(z,  q)  =  p(z)  +  qa(z)  be  the  characteristic 
polynomial  of  an  explicit  LMM  where  p  and  a  are  of  degree  k.  and  k-1, 
respectively.  The  polynomial  1/q  n(z,  q)  has  the  same  roots  as 
ll(z,  q).  We  need  only  consider  real  q  and  the  coefficient  of  z 


in  1/q  !!( 2 ,  q)  us  q  approaches  negative  infinity.  The  previous 
lemma  guarantees  an  N  such  that  ii(z,  q)  has  roots  exceeding  unit 
magnitude  for  q  <  N.  QLD. 

Jeltsch  [20]  finds  that  Cryer's  [4]  methods  are  not  only  A^-- 
stable  but  also  A(0)-stable.  Hence  there  exists  A(Q)-st«b'!e  methods  cl 
arbitrarily  nigh  order.  Gupta  [14]  has  found  A(  a) -stable  me  i  ho  Jr. 
through  order  12  with  a  exceeding  70° ,  however,  the  truncation  on c-r 
gets  extremely  large  on  the  higher  order  methods. 

Another  useful  definition  is  provided  by  Oc-cn  and  Linger  [30]. 

Definition:  A  LKM  is  A  -stable  if  'it  is  absolutely  stable  in 
a  neighborhood  of  infinity  on  the  complex  q- soil  era.  We  take  the 
following  theorem  from  Harden  [29]  without  proof. 

Theorem  2.3-4  (Rouchd).  If  P(z)  and  Q(z)  are  analytic  interior 
to  a  simple  closed  Jordan  curve  C  and  if  they  arc  continuous  on  C  and 


|P(z)|  <  |Q(z)j  for  zcC,  then  the  function  F(z)  =  P(z)  +  Q(z)  has  the 
same  number  of  zeros  interior  to  C  as  does  Q(z). 

The  following  lemma  is  based  on  Rouchd's  Theorem  and  is  ussfel 
in  providing  insight  into  the  concept  of  A  -stability.  In  the 

remainder  of  this  section  let  L  be  an  implicit  zero-stable  L MM  with 

characteristic  polynomial  h'(z,  q)  =  p(z)  +  qa(z)  and  C  =  { z j  J z |  =  1). 

Lemma  2.3-5.  Consider  L  as  described  above  and  suppose  none 
of  the  roots  of  c  lie  on  C.  Then  there  exists  a  real  number  S  such 

that  for  | q |  >  S,  li(z,  q)  and  o(z)  have  the  same  number  of  zeros 

interior  to  C. 


Proof:  Let 
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S  = 


sup 

zeC 


|p(  z) 

lotzy 


Then  S  is  finite  and  if  |qj  >  S  we  have  |qa(z)|  >  |p(z)|  for  zeC. 

By  Rouchd's  Theorem  we  know  that  3l(z,  q)  and  qa(z)  have  the  same 
number  of  zeros  interior  to  C.  Since  the  roots  of  qa(z)  and  o(z) 
are  identical  and  since  the  degrees  of  II  and  a  are  the  same,  we  have 
Il{z,  q)  and  a( z)  with  the  same  number  of  zeros  interior  to  C.  QED. 

Now  we  prove  a  simple  and  useful  characterization  of  A^-stable 

methods. 

Theorem  2.3-6.  L  is  A  -stable  if  and  only  if  all  the  roots 
of  a(z)  lie  inside  C. 

Proof:  Let  the  degrees  of  n  be  k  and  let  S..  and  z^q), 
i  =  1 ,  ••  •  ,  k  be  the  roots  of  o(z)  and  Il(z,  q),  respectively, 
where 


lim  zi (q)  =  S 

q-K» 

If  L  is  A^-stable  there  is  a  neighborhood  of  q  =  «  in  which 
{ (q ) |  <  1  so  that  clearly  j$. |  <  1 ,  i  =  1 ,  •••  ,  k. 

If  |S.j  |  <  1 ,  i  =  1 ,  * •*  ,  K  then  the  previous  lemma  provides 
A^-stabi lity  immediately.  QED. 

The  result  which  follows  provides  a  relationship  between 
Ag-stabilizy  and  A^-stabil  ity. 

Theorem  2.3-7.  If  L  is  Aq-s table  and  none  of  the  roots  of 
o{z)  lie  on  C,  then  L  is  A  -stable. 
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Proof:  Suppose  one  of  the  roots  of  a,  say  S .  has  magnitude 

J 

greater  than  one.  The  previous  lemma  provides  a  real  number  S  with 
a  root  of  Jl(z,  q^)  outside  C  for  all  q^  with  | |  >  S.  We  therefore 
contradict  A^-stability  and  we  conclude  all  roots  of  a  lie  inside  C 
and  hence  L  is  A  -stable.  QED. 

The  last  definition  we  look  at  from  the  literature  is  from 
Gupta  [13]  and  combines  features  of  Gear's  stiff  stability  and  the 
A(a)-stability  of  Widlund. 

Defi ni tion:  A  LMM  is  said  to  be  A(ct,  D)-stable,  ctc(0,  TI/2) 
if  the  region  of  absolute  stability  includes  all  q  with  |Arg(-q)|  <  a, 
q  ^  0  and  all  q  with  Re  q  _<  D. 

Of  all  the  definitions  reviewed  above,  the  only  one  to  define 
a  region  of  accuracy  about  the  origin  was  Gear's  stiff  stability 
definition.  As  mentioned  before  the  reference  to  accuracy  in  that 
definition  is  vague.  This  has  led  to  different  interpretations  of 
what  is  meant  by  the  region  of  accuracy  about  the  origin  for  example 
see  Oeltsch  [20,  p.  9]  and  Gupta  [15,  p.  492].  The  motivation  for 
developing  stiffly  stable  methods  was  to  allow  us  _o  determine 
stepsize  on  the  basis  of  accuracy  constraints  alone.  To  that  end 
we  do  need  a  description  of  the  infinite  region  of  stability  as 
provided  by  the  foregoing  definitions.  However,  we  cannot  overlook 
then  the  parallel  need  to  provide  a  description  of  the  region  of 
accuracy  with  equivalent  precision.  Absolute  stability  as  mentioned 
before  provides  a  measure  of  accuracy  which  is  too  lax  in  the  left 
half  plane  and  too  rigid  in  the  right  half  plane.  We  must  be 
concerned  with  digits  of  precision  when  performing  computations  on  a 
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digital  computer.  Actually  digits  of  precision  are  a  measure  of 
relative  error.  For  these  reasons  then  it  seems  natural  to  continue 
our  concern  about  relative  error  with  a  measure  of  relative  stability 
which  we  include  in  the  definition  given  below. 

Definit.i on :  A  LMM  is  said  to  be  A(a,  r)-stable  if  the 
method  is  A(a)-stable  with  regard  to  its  region  of  absolute  stability 
and  is  relatively  stable  within  the  disk  of  radius  r  about  the 
origin. 

It  is  this  definition  which  v/e  implement  in  our  search  for 
stiffly  stable  LMM.  A  method  which  is  A(a,  r)-stable  will  be 
relatively  robust  depending  upon  the  size  of  a  and  r.  By  that  we 
mean  it  should  give  good  results  on  a  large  variety  of  probleins--both 
stiff  and  non- stiff. 


CHAPTER  III 


ERROR  ANALYSIS 

3.1  Introduction 

Error  in  our  approximation  to  the  solution  of  (1.1-1)  occurs 
because  we  base  our  recursive  calculations  on  inaccurate  values, 
because  our  recursive  model  does  not  accurately  represent  the  problem, 
and  because  our  calculations  themselves  are  not  always  precise.  Also, 
the  recursive  nature  of  our  model  inherently  propagates  existing  error. 
Controlling  this  propagation  of  error  is  the  aim  of  stability 
constraints  discussed  in  Chapter  II.  Any  complete  and  careful 
analysis  of  cumulative  error  must  consider  at  least  these  sources  of 
error. 

3.2  A  Global  Sound 

Henri ci  [16,  section  5.3-4]  investigates  this  cumulative 
error  and  arrives  at  the  bound  given  below  for  the  solution  of 
(1.1-1)  by  a  member  of  C(p,  k).  If 

|  of  j  6_i |hL  <  1 

we  have  the  global  error  en  satisfying 

ienl  1  T*[A6k  +  (xn-a)(k^  +  GyhP)]exp[(xn-a)  LT*B]  (3.1-1) 
where  the  following  definitions  are  used. 


k-1  k-1 

A*  I  !«,!.  B  =  I  13-1, 

i=-l  1  i=-l  1 

L  is  the  Lipschitz  constant  for  the  function  f,  6  is  the  maximal 
starting  value  error, 

r  =  sup  (lAjD^o 

where 


1  r  ,  i 

m =  1 


and 


k-1 

Hz)  =  l  a  .2 

i=-l  1 


i+1 


(r  is  shown  to  be  finite  in  Henri ci  [16,  P-  242]  for  zero-stable 
methods. ) 


r*  = 


l-h|a_-j3  ^|L 

h  is  the  stepsize,  xne[a,  b],  is  a  bound  for  the  magnitude 

of  the  roundoff  error  committed  at  any  step  of  the  integration  process, 


G  =  /  |G(s) 1 ds 
0 
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where 


G(s 


k-1 

)  =  -Li  ak-j-2^J'_s^+  “  P0k-j-2^J“S^ 

3  **  • 


and 


(j-s)+  =  max  {0,  j-s), 
and 

y  =  max  |t/p+1  ^  (x)  j . 
xe[a,b] 

Error  bounds  generally  have  the  burden  of  coping  with  worst  case 
examples  and  therefore  many  times  are  forced  to  give  outrageously 
large  bounds  on  simple  problems.  This  severely  damages  their 
practical  value. 

Consider  for  example  the  bound  (3.1-1)  applied  to  the 
solution  of 


{/'  =  A<y,  y( 0)  =  1,  xe [0,  1]  (3.1-2) 

using  the  second  order  backwards  difference  formula 

Vi  -  «/3  *.  -  >/3  Vi  +  2/3  hf'Vi-  W  <3J-3> 


p(z)  =  -  z2  +  4/3  z  -  1/3 


with  stepsize  h.  Then 


so  that 


1 

Jft T 


For  ] z |  <  1  we  find  the  Laurent  expansion 


1  .  v  [l-3Wlj 

£(z)  1=0  2-31 


Hence 


n  **i+l 

T  =  sup  { 1  A.  | }°°  =  sup  {  ■— --•■■  }?_  n  =  3/2. 

1  i=0  '  2*3  1  u 


In  general  for  the  method  (3.1-3)  we  see  that  (3.1-1)  reduces  to 


i  n  ?  DL(x • 

K1  -P3hl  M85  +  <Va><9k/  +  2yh£)J  exp  3.^L'-  |  (3.1-4) 


and  is  valid  as  long  as  hL  <  3/2. 


In  this  example  we  then  have  the  bound  as 


485  +  9n  k,hq+1  +  2nLV 
!enl  ± - 6^4hL -  exP 


3nhL  { 
3-2hy ' 


(3.1-5) 


If  A  =  -100  and  h  =  0.01  then  hL  =  1  <  3/2  so  that  this  expression  is 


valid.  We  compute  this  bound  on  the  error  after  25  steps  and  in 


doing  so  it  is  of  negligible  consequence  to  omit  the  positive  contri¬ 


butions  of  starting  error  and  roundoff  error.  In  this  case 


me 
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I I  <  25  exp  (75)  £  9.33  •  10^  when  the  true  solution  {/(Xgg)  ^ 

.138  •  TO"10.  Clearly  the  information  provided  by  this  bound  from  a 
practical  point  of  view  is  totally  useless.  When  error  bounds  behave 
in  this  way  we  can  sometimes  get  better  results  from  the  use  of  error 
estimates  [9,  p.  14],  However,  as  with  the  example  given  by  Gear  [9, 
p.  16],  error  estimates  can  give  smaller  numbers  than  the  actual  error. 

From  a  practical  point  of  view  it  would  be  more  useful  if 
the  expressions  for  roundoff  error  and  starting  error  were 
respectively  dependent  upon  machine  precision  and  accuracy  of  the 
starting  procedure.  The  bound  (3.1-1)  is  theoretically  satisfying 
from  the  viewpoint  that  it  does  approach  zero  with  h.  It  would  be 
pleasing  to  see  the  effect  of  stability  properties  represented  in 
the  bound.  In  the  next  section  we  consider  a  representation  of  the 
global  error  which  is  more  useful  for  certain  practical  and 
theoretical  considerations. 

3.3  An  Approximate  Error  Representation 

As  in  section  1.2,  if  we  assume  our  solution  {/(x)  has 
sufficient  derivatives,  the  local  truncation  error  can  be  expressed 
as 


L(K(xnt1),  h)  •  <W,  + 

-  Cp+,  h1*1  !/(P+,!(x„<.1)  +  0(hp<'2)  (3.3-1) 


for  a  member  of  C(p,  k)  where  this  notation  is  taken  from  Henri  ci 
[16,  p.  220].  L,  as  defined  above,  may  be  regarded  as  a  linear 
operator  on  any  differentiable  function  i/(x),  and  if  y  possesses 
sufficient  derivatives  its  value  can  be  given  as  the  right  side  of 
(3.3-1).  Following  the  investigation  of  Henrici  for  arbitrary 
operators  of  the  form  (3.3-1),  others  [27,  33]  have  pointed  out  that 
in  all  cases 


|L(i)(xnt1).  h)l  <  ^rhptlGy,  (3.3-2) 

where  G  and  y  are  defined  as  in  the  previous  section.  G(s),  also 
defined  in  the  previous  section,  is  called  the  influence  function. 

In  those  cases  where  G(s)  is  of  the  same  sign  throughout  the  interval 
[-1,  k-1],  we  have 

lMxnt)),  h)  =  Cp+1  h^1  6r(p+1)(C).  (3.3-3) 

for  €e(x  _k+-|,  xn+])*  G(s)  does  not  change  sign  on  [-1,  k-1]  for 

many  commonly  used  methods  including  Adam's  formulae  and  the  backwards 

difference  formulae  of  Gear  with  order  below  six. 

We  now  develop  an  approximate  representation  of  the  global 

error  which  occurs  in  the  solution  of  yl  =  hj  by  a  member  of  C(p,  k). 

Let  i/(xn)  and  yn  be  the  true  and  approximate  solution,  respectively. 

Let  the  local  truncation  error  and  roundoff  error  occurring  in  the 

calculation  of  y  be  T  and  R  respectively.  Then  T  =  L(t/(x  ),  h) 

■'n  n  nr  J  n  -  n 


k-1 


t 

■>  _  V 

V.  "  .L 


W,  *  h-v,-  )/,._  i  ,  • 


i  -  r:- 


If  we  define  c  =  :/{x  }  -  y  wc  have 

fi  "  fi  Ft 


T  .  -,  -  K  i 
n-<-l  (!■>  1 


fc-1 

r 


r-. 


hx?,.  }'.{?■  - 

;  n- : 


k- 1 

\  4  ■  *  ^  "»  \ 

l  ^  -  n -k . /y  . 

-»  I  «  « !-  ! 

i - I 


x-1 

I 

• _  •% 

1  —  s 


(2-  -s  h; '  )*2  . 

i  1  n-i 


{3. 3-4  * 


If,  as  in  Lor.terl  [27,  0.  651,  we  assume  T  0  - 

n  "  1  n  •  i 

we  have 


a  C iSfil  > 


.  _  c  J  .n  ' 
e  =  )  d.Z .  t  -. ** ; — 

n  1  1  i.'.ot-- 


{3.3-5:. 


1=0 


where  z.  ere  the  roots  of  .he  characteristic  polyiionnai  11(2.  ht). 
This  results  as  the  solution  to  the  non- horr.oqenecus  constant 


When  solving  a  problem  such  us  (l.l-i)  with  values  of  hi 

inside  of  absolute  stability  but  outside  the  region  of  relative 

stability,  the  miricatieu  or* or  will  approach  zero  but  rot  ss  rapidly 

as  the  true  solution.  The  calculation  then  is  fenced  to  ccrcim.;- 

un til  the  error  coa'os  within  accuracy  constraints  even  though  the 

actual  solution  nvy  ha'-.,  long  been  satisfactorily  small.  The  expense 

of  this  continued  calculation  offsets  and  conceivably  could  exceed 

the  expanse  of  using  a  smaller  slepsize  to  Dring  hA  within  the  region 

of  relative  stability.  This  is  why  relative  stability  can  be 

important  even  in  the  left  half  plane.  Examples  of  this  were 

obtained  in  computer  runs  using  3b- dig  it  precision.  Solution  of 

(3.1-2)  with  A  -  100  and  A  =  -1  using  the  method  (3  1-3)  and  h  ::  0.01 

8  -  5 

gave  relative  error  of  2.33  •  10  and  3.3  •  10  ,  respectively,  after 

100  steps. 

Inside  the  region  of  relative  stability  the  second  term  in 
(3.3-5)  usually  dominates  the  global  error.  The  problem  independent 
part  of  the  local  and  global  truncation  error  constants  are 


3*1  and  tCl 

X  -  0(1)’ 


respectively.  If  we  use  sufficient  machine  precision  to  insure  the 
roundoff  error  is  negligible  compared  to  the  truncation  error,  then 
in  (3.3-5)  <f  is  dominated  by  the  local  truncation  error 


Cp+1 


Tills  robes  it:  appaivn!"  why  the  pivSic!:*  ! ude, 'end--::!  pan-  of  the  local 
truncation  error  diffr'-t  from  the  global  ininca- io"  error  by  a 
factor  of 


This  observation  is  supported  by  computer  runs,  ay., in  in  35- digit 
precision,  where  in  the  sciatic::  of  (3.1-f.)  with  A  -  -1  using  a 
variety  of  LKii  and  stopsir.es  the  latio  cf  local  to  global  trur.caiion 
error  was  consistently 


within  approximately  three-  significant  digits.  It  is  important  to 
note  here  that  there  is  an  entire  subclass  of  C(r,  k)  which  have  local 
truncation  error  constants  larger  than  their  global  truncation  error 
constants,  and  conversely.  This  results  since  the  determination  of 
Cp+i  and  o(l)  are  independent  in  the  solution  of  the  system  of 
equations  (1.2-3)  which  defines  the  members  of  C(p,  k).  Evaluation 
of  the  relative  merits  of  different  LMM  need  to  take  this  into  account. 
In  the  investigation  reported  in  Chapter  V  the  term 


of  IT 


is  used  as  the  measure  of  accuracy  for  comparison  of  LMM. 


CHAPTER  IV 


PREDICTOR  SEARCH 


4 . 1  Statement  of  and  Approach  to  the  Problem 

Predictor  methods  commonly  suffer  smaller  stability  regions 
and  larger  truncation  error  constants  than  implicit  nethods  [27, 
p.  84].  On  the  other  hand,  implicit  methods  suffer  the  difficulty 
of  finding  the  solution  to  the  implicit  relation  defining  the  forward 
solution  approximation  at  each  step  of  the  integration.  The 
solution  to  these  problems  is  normally  accomplished  by  pairing  a 
predictor  with  a  corrector  to  form  a  predictor-corrector  ( P-C) 
algorithm.  In  the  solution  of  non-stiff  problems,  if 
{yi  |  i  =  0,  ,  N)  is  our  discrete  model  of  the  solution  to 

(1.1-1),  the  predictor  is  used  to  compute  a  first  approximation  to 
forward  values  yp+-j-  This  predicted  value  is  then  used  as  yjj]:j  in 
(1.2-5)  to  start  the  recursive  evaluation  of  the  sequence 


(vtsV 

lyn+Vs=l 


This  sequence  may  be  evaluated  as  far  out  as  desired  although 
convergence  within  desired  accuracy  constraints  is  usually  very 
fast.  For  example  the  popular  code  DIFSUB  by  Gear  [10]  restricts 
s  _<  4.  In  the  solution  of  stiff  problems,  the  condition  required 
for  convergence  of  the  corrector  iteration. 
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IhB^ofJlL  <  1 

where  L  is  a  Lipschitz  constant  for  f(x>  tj)  with  respect  to  t/,  forces 
us  to  discard  corrector  iteration  in  fa-'or  of  a  Newton  iteration 
[27,  p.  239].  When  Newton's  method  is  applied  to  (1.2-5)  for  the 
solution  yn+i ,  the  restriction  on  h  required  for  convergence  is  usually 
much  more  relaxed  than  that  imposed  by 

|h0_^a_]|L  <  1. 

But  even  in  stiff  problems,  predictors  are  used  to  provide  accurate 
initial  estimates  for  the  Newton  iteration. 

If  predictors  with  significantly  improved  stability 
characteristics  exist,  at  least  two  interesting  possibilities  arise. 
Those  possibilities  are  improved  P-C  algorithms  and  cedes  using 
predictors  exclusively.  In  a  P-C  algorithm  as  described  in  the 
previous  paragraph,  unless  the  corrector  iteration  is  carried  out  to 
convergence,  the  stability  characteristics  of  the  P-C  algorithm  are 
not  those  of  the  corrector  alone  [1].  The  fewer  the  number  of 
corrector  iterations  performed,  the  less  heavily  the  s tabi lit"  of 
the  P-C  algorithm  depends  on  the  corrector.  Since,  as  noted  above, 
in  practice  few  iterations  are  actually  carried  out  it  leads  us  to 
believe  that  it  may  be  important  to  use  a  predictor  with  as  large  a 
stability  region  as  possible. 

The  other  possibility  of  using  predictors  alone  for  non-stiff 
problems  is  very  attractive  in  all  cases  and  occasionally  the  use  of 


explicit  methods  is  necessary.  Among  the  advantages  to  be  gained  are 
economy,  simplicity,  and  versatility.  Economy  is  realized  with 
regard  to  speed  of  computation,  programming  efforts  in  developing 
a  reliable  code,  and  machine  storage.  For  most  non-trivial  problems, 
by  far  the  most  expensive  computation  performed  in  finding  the  model 
is  that  of  function  evaluations  for  the  derivative  approxinia- 
tions  y'+.].  Predictor  algorithms  reduce  this  to  one  function  evaluation 
per  step  whereas  P-C  algorithms  require  from  two  to  four.  Elimina¬ 
tion  of  the  corrector  iteration  and  the  associated  coding  greatly 
reduces  the  programming  effort.  It  also  simplifies  considerations 
such  as  stability  analysis.  For  example  see  Lambert  [27,  p.  97]  for 
a  definition  of  the  characteristic  polynomial  of  a  P-C  method.  A 
proper  stability  analysis  must  examine  the  roots  of  the  characteristic 
polynomial  which  we  see  is  much  more  complicated  than  for  a  predictor 
alone.  The  reduced  coding  and  storage  required  of  predictor  algorithms 
render  such  codes'  implementation  within  the  capabilities  of  smaller 
machines,  even  hand-held  calculators  when  the  problem  is  not  large. 

Gear  [11]  gives  real-time  integration  as  an  example  of  where 
implicit  methods  cannot  be  used  in  the  usual  sense.  Here  again  the 
use  of  predictors  would  be  natural  were  it  not  for  the  smaller  stability 
regions  of  these  explicit  methods. 

It  is  true  that  normally  predictor  methods  also  have  larger 
truncation  error  constants  than  correctors  but  this  can  always  be 
more  than  compensated  for  by  using  reduced  steplength.  For  example 
in  comparing  the  truncation  error  constants  of  the  Adams -Bash forth 
to  Adams-lloulton  methods,  the  stepsize  reduction  factor  required  to 
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equilibrate  the  truncation  error  constants  for  a  given  order 
increases  from  0.447  for  second  order  to  0.671  for  ninth  order.  It 
is  the  smaller  stability  regions  which  create  the  greatest  drawback 
for  the  exclusive  use  of  predictors  in  solving  certain  problems 
(1.1-1). 

To  date  no  comprehensive  search  has  been  conducted  for  more 
stable  predictors.  Those  searches  existing  in  the  literature  have  had 
limited  specific  goals,  for  example  see  [24,  3].  The  purpose  of  the 
investigation  reported  in  this  chapter  then  is  to  determine  how  far 
the  stability  regions  of  predictor  methods  can  be  extended.  These 
methods  may  find  utility  by  incorporation  into  P-C  algorithms  or  in  a 
code  which  employs  explicit  methods  exclusively.  (A  related  question, 
not  investigated  here,  is  whether  or  not  it  is  the  case  that  P-C 
algorithms  with  better  stability  properties  are  always  the  result  of 
combining  separate  predictors  and  correctors,  each  with  good  stability 
properties. ) 

A  guide  as  to  the  class  of  predictors  in  which  we  may  expect 
favorable  results  is  provided  by  Henrici  [16,  section  5.2-8].  The 
theorem  is  given  here  without  proof. 

Theorem  4.1-1 .  The  order  of  a  zero-stable  LMM  whose 
stepnumber  is  k  cannot  exceed  k+1  when  k  is  odd  or  k+2  when  k  is  even. 

As  discussed  in  Lambert  [27,  p.  57],  those  zero-stable 
methods  of  the  class  P(k+2,  k)  when  k  is  even  have  no  interval  of 
absolute  stability.  We  ere  thus  restricted  to  investigation  of 
classes  P(n,  m)  where  n  <  m  +  1.  We  choose  the  class  P(n,  n)  because 
of  the  additional  degree  of  freedom  permitted  over  P(n,  n-1)  and 
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because  many  codes  already  in  existence  have  implemented  the  Adams- 
13  ash  forth  methods  which  lie  in  this  class. 

We  use  the  relative  stability  region  of  a  LMM  as  our  measure 
of  stability.  We  could  use  absolute  stability  regions  for  such  a 
measure,  however  this  would  provide  us  no  knowledge  of  the  behavior 
of  our  approximate  solution  to  (1.1-1)  for  q  in  a  neighborhood  of 
the  origin.  We  define  the  radius  of  relative  stability  to  be  the 
radius  of  the  largest  circle  centered  at  ti.e  origin  and  contained  in 
the  region  of  relative  stability. 

There  are  several  definitions  of  relative  stability  existing 
in  the  literature  [27,  p.  68],  and  it  is  important  to  distinguish 
which  definition  is  used  when  making  comparisons  of  stability  regions 
from  different  sources.  Even  though  for  small  q  we  have  the 
relationship  (2.1-1)  there  can  be  large  differences  for  q  near  zero. 

In  our  investigation  we  use  the  definition  given  in  section  2.2  which 
is  taken  from  Lan&ert  [27].  Another  approach  to  defining  relative 
stability  is  to  compare  the  magnitude  of  the  roots  of  the  charac¬ 
teristic  polynomial  to  the  magnitude  of  exp  (q)  rather  than  to  the 
magnitude  of  the  principal  root.  For  example,  Crane  and  Klopfenstein 
require  | rg  |  _<  exp  (q) ,  s  =  1 ,  2,  •••  ,  k  and  that  roots  of  magnitude 
exp  (q)  be  simple.  For  comparison,  consider  the  second  order  backwards 
difference  formula  (3.1-3).  This  method  has  characteristic  polynomial 
n(z,  q)  =  2/3{q-l)z^  +  4/3  z  -  1/3.  For  the  sake  of  simplicity,  we 
restrict  our  comparison  to  the  interval  of  relative  stability  which 
is  the  real  values  of  q  for  which  the  method  is  relatively  stable. 
Lambert's  and  Crane  and  Klopfenstein's  definition  yield  the  values 


(-1/2,  3/2)  U  (3/2,  ®)  and  (a,  0],  respectively,  where  a  %  -0.75262. 

Other  methods  tested  indicate  possibly  for  Crane  and  Klopfenstein's 

definition  the  radius  of  the  complex  region  of  relative  stability  is 

always  zero.  For  example  this  is  the  case  when  the  boundary  of  the 

region  of  absolute  stability  intersect.-  imaginary  axis  only  at  the 

origin.  This  follows  since  for  q  imagina.  ,  {exp  (q)|  =  1,  and  the 

condition  of  relative  and  absolute  stability  coincide.  Also  if  ij  and 

y  are  the  true  and  computed  solutions  at  a  point  x  ,  then  the  latter 
n  n 

definition  restrains  jyj  =  \ij  |  -  c,  for  c  >  0.  The  magnitude  of  c 
and  not  its  sign  determine  the  acceptability  of  y^.  For  these 
reasons  we  find  the  latter  definition  unsuitable.  Lambert's  definition 
ties  the  growth  of  the  relative  error  to  the  size  of  the  principal 
root,  however,  for  q  not  large  and  especially  for  higher  order 
methods. (2. 1-1 )  indicates  this  is  not  a  lax  requirement.  Sons 
definitions  of  relative  stability  which  compare  magnitudes  of  the 
roots  of  the  characteristic  polynomial  to  that  of  the  principal  root 
do  not  require  strict  inequality,  see  for  example  Ralston  [33,  p.  177]. 
However,  strict  inequality  more  easily  lends  itself  to  computer 
implementation. 

We  need  then  to  find  a  means  of  evaluating  the  radius  of 
relative  stability  defined  above.  An  expression  for  the  boundary  of 
the  region  of  relative  stability  cannot  be  found  as  with  absolute 
stability.  It  is  necessary  to  actually  find  the  roots  of  the 
characteristic  polynomial.  We  define  a  function  introduced  by 
Rodabaugh  [34]  and  called  the  critical  difference  function.  Given 
a  LMMwith  characteristic  polynomial  Jl(z,  q) ,  we  define 
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C(q)  =  |z1  (q) |  -  max  { 1 2$(q) i >s=2 

where  (q) ,  S  =  1,  2,  •••  ,  k  are  the  roots  of  JI(z,  q)  and  z-j(q) 
is  the  principal  root.  In  our  search  for  stable  predictors  we 
require  C(0)  >  0,  or  strict  zero-stability.  We  are  thus  guaranteed 
a  neighborhood  of  stability  about  the  origin  by  the  continuity  of 
the  critical  difference  function. 

To  determine  the  radius  of  relative  stability  we  need  then 
to  locate  the  zero  of  C(q)  nearest  the  origin.  The  radius  will  be 
equal  to  the  modulus  of  this  zero. 

We  do  this  numerically  by  taking  a  finite  set  of  rays  through 
the  origin,  finding  the  zero  of  C(q)  nearest  the  origin  on  each  ray 
in  this  set,  and  approximating  the  radius  as  the  smallest  modulus 
possessed  by  the  resulting  zeros.  However,  to  locate  the  zero  of  C(q) 
nearest  the  origin  on  a  given  ray  we  must  be  able  to  evaluate  C(q) 
and  this  entails  tracking  the  principal  root. 

To  track  the  principal  root  along  a  ray,  q  is  given  a  fixed 

argument  and  its  modulus  is  incremented  successively  by  a  small 

amount.  We  know  z,(0)  =  1  and  if  q  and  q  are  successively 
I  ti  w  1 

incremented  values  of  q  on  a  given  ray,  then  z^ (qn+-j )  is  assigned  to 
be  the  root  of  H(z^  qn+-j)  nearest  z^(qn).  Thus  starting  at  the 
origin  C(q)  is  evaluated  at  the  successive  increments  of  q  as  we 
track  the  principal  root.  The  first  increment  at  which  C(q)  is 
negative  identifies  an  interval  which  contains  the  first  zero  of 
C(q).  A  bisection  routine  is  then  used  to  more  accurately  identify 
the  zero.  It  is  conceivable  that  the  first  zero  of  C(q)  could  be 
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missed  in  this  way,  however,  it  would  be  necessary  for  it  to  be 
missed  on  several  rays  in  order  for  the  output  value  of  the  radius 
to  be  severely  affected.  This  is  highly  unlikely  but  we  should  remain 
aware  of  the  possibility. 

Evaluation  of  the  radius  of  relative  stability  can  be  costly 
and  rather  complex.  It  is  possible  that  much  research  into  relative 
stability  characteristics  possessed  by  LKM  has  been  inhibited  by 
these  factors.  In  addition,  these  factors  may  have  partially 
moti rated  the  alternative  definitions  which  compar.  the  roots  of 
n(z,  q)  to  exp  (q)  rather  than  tc  the  principal  root.  There  are 
several  considerations  which  greatly  reduce  the  cost  factor  involved. 

These  considerations  are  based  on  properties  of  the 
characteristic  polynomial  fi(2,  q)  and  the  resulting  boon da  17  of 
consistent  zero-stable  methods.  First  of  all,  it  suffices  to  check 
for  zeros  of  C(q)  when  Iin  (q)  _>  0,  since  Lemma  2.2-1  provides  that  the 
boundary  of  the  region  of  relative  stability  is  symmetric  about  the 
real  axis.  We  also  economize  by  restricting  Re  (q)  _<  0  since  in  practice 
it  is  there  that  we  find  the  most  severe  restrictions  on  the  radius. 
Normally  we  find  it  sufficient  to  determine  zeros  of  the  function  C(q) 
on  15  rays  in  the  second  quadrant,  however,  occasionally  for  verifi¬ 
cation  we  increase  this  number  and  extend  our' investigation  "into  the 
first  quadrant.  Since  we  search  only  for  the  radius,  we  terminate 
our  search  for  a  zero  of  C(q)  on  a  given  ray  once  we  have  incremented 
the  modulus  of  q  bevond  the  modulus  of  any  previously  determined 
zero.  Thus  due  to  the  shape  of  a  typical  boundary  it  is  less  time 
consuming  to  assign  the  first  ray  an  argument  of  n  and  proceed 


backwards  to  an  argument  of  H/2  rather  than  conversely.  Finally,  based 

on  (2.1-1)  we  expect  exp  (q)  =  z^(q)  for  q  small.  A  computer  analysis 

of  this  approximation  has  shown  this  to  be  of  practical  value  as  well. 

For  a  large  number  of  methods  tested  through  ninth  order,  when 

| q |  £  0.25  the  different  jexp  (c)  -  2^ (q ) |  was  always  at  most  on  the 
_2 

order  of  10  .  In  practice  then  we  use  the  alternate  means  of 
identifying  the  principal  root  as  that  root  of  K(z,  q)  nearest  exp  (q), 
for  ( q 1  <  0.25.  In  this  range  of  moduli  much  larger  increments  are 
permitted.  This  feature  alone  lias  resulted  in  a  cost-savings  of 
approximately  The  savings  here  is  greater  when  the  radius  of  the 

method  is  small,  for  example  in  higher  orders.  It  is  difficult  to 
assess  the  reduction  in  the  cost  of  this  investigation  afforded  by  trie 
combination  of  factors  mentioned  above,  however,  it  is  clear  that  they 
have  vastly  increased  its  feasibility  from  the  cost  perspective. 

The  search  for  stable  predictors  within  the  class  P(n,  n) 
is  posed  as  an  optimization  problem.  In  order  to  do  this  we  isolate  a 
specific  real  quantity  which  is  to  be  minimized.  That  quantity  is 
the  negative  of  the  radius  of  relative  stability.  Clearly  there  are 
limitations  incurred  by  selecting  a  specific  objective  such  as  we  have. 
For  example  in  the  nature  of  the  resulting  coefficients  and  error 
constants  of  the  optimal  method.  However,  our  results  prove  valid  as 
inputs  to  decision  processes  concerning  further  development  of  numerical 
methods  for  approximating  the  solution  of  problems  (1.1-1),  and  as 
numerical  methods  in  their  ov/n  right  for  some  classes  of  problems. 
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Two  packaged  codes  were  incorporated  into  the  program 
written  to  locate  stable  predictors.  The  optimization  code  used 
(ZCXPWLD)  is  based  on  Fletchers  [8]  version  of  Powell's  [32]  penalty 
function.  It  uses  Powell's  unconstrained  optimization  algorithm  [31] 
to  determine  a  local  minimum  through  a  sequence  of  unconstrained 
minimization  problems.  For  the  most  part  it  was  not  necessary  to 
use  constraints. 

In  those  cases  we  gained  convergence  more  rapidly  and  the 
algorithm  reduced  to  that  described  in  [31].  .Another  code  (ZCPOLY) 
employed  in  our  search  finds  the  zeros  of  a  complex  polynomial.  This 
code  is  in  the  International  Mathematical  and  Statistical  Libraries 
(IMSL).  It  uses  the  Jenkins-Traub  three  stage  variable  shift 
iteration  described  in  [21]  and  L22].  The  remainder  of  the  program 
was  developed  by  this  writer. 

The  domain  of  our  object  function  for  the  class  P(n,  n)  is 
pP~\  This  results  from  consideration  of  the  system  of  equations 
(1.2-3)  which  determines  the  class  P(n,  n)  as  discussed  in  section 
1.3.  The  only  modification,  to  include  notation,  of  that  discussion 
applicable  here  is  that  for  predictors  e  ^  =  0.  The  flow  of  the  program 
consists  of  taking  a  point  (a^ ,  •••  .  an  in  Rn-^  and  determining 
the  coefficients  of  the  corresponding  LMM  according  to  the  procedures 
outlined  in  section  1.3.  We  then  evaluate  C(0).  If  C( 0)  _<  0  the 
object  function  is  assigned  a  zero  value.  Otherwise  the  radius  is 
evaluated  as  previously  described  and  the  object  function  is  assigned 
the  negative  of  the  radius.  Based  on  this  assigned  function  value. 


1 

the  optimization  algorithm  calculates  a  new  domain  point  in  k  . 

We  then  repeat  this  cycle  until  functional  values  converge  within  a 
prescribed  tolerance.  This  optimization  code  perform  '  veil,  usually 
converging  within  100  function  evaluations. 

Essential  for  convergence  is  a  starting  point  (method)  with 
a  positive  radius.  Good  starting  points  are  readily  available  from 
the  literature,  for  example  the  Adams-Bashforth  methods,  lie  may  gain 
additional  starting  values  from  consideration  of  a  theorem  from 
Chapter  II.  K  we  factor  p(z)  =  (z-l)p^(z),  then  those  polynomials 
p(z)  for  which  p^(z)  is  a  Schur  polynomial  are  associated  with  LKM  of 
positive  radius.  Theorem  2.1-4  provides  a  characterization  of  the 
smallest  convex  set  in  H‘  containing  the  coefficients  of  all  Schur 
polynomials  of  degree  n.  For  example  we  compute  this  characterization 
explicitly  fcr  information  regarding  starting  value:  for  the  class 
P(5,  5).  We  have  a  -1  and 


“o  “ 1  -  .1,  °i 


so  that  p(z)  =  (z-1)  p^(z),  where 


Ps(z)  =  I 

*  i=0 


i=0 


4-i 


We  apply  Theorem  2.1-4  to  the  polynomial  p~(z). 


The  matrix  is  given  by 


1-1  1-11 
-4  2  0  -2  4 


0  -2  0  6 


(a0,  ar  a2,  a3,  a4)  = 


-4  -2 


1111 


(4.1-1) 


(«-,»  a3  +  «2  +  '3  "  a4’  al  +  a2  +  a3  +  a4’  1) 


The  necessary  condition  in  the  Theorem  is  satisfied  by  requiring  all 
elements  of  the  product  M4(aQ,  ***  ,  an) 1  to  be  positive.  This  leaves 
us  with  the  following  systems  of  inequalities. 


m 


1  -  a-j  -  «3  >  0 

2  -  a.  -  a„  -  2a.  >  0 

12  4 

3  “  c2  “  a3  +  ^a4  >  ® 

2  +  a.  +  a0  -  2a,  >  0 

12  4 

1  +  Oj  +  20^  +  3a^  +  4a^  >  Q. 


(4-1-2) 


4.2  Results  and  Comparisons 

The  results  from  the  investigation  described  in  the  previous 
section  are  given  here  for  orders  four  through  nine.  Since  the  Adams - 
Bashforth  (A-B)  methods  are  used  nearly  exclusively  in  popular  codes 
employing  explicit  UK,  and  since  the  K-th  order  A-B  method  is  a 
meiriber  of  P(K,  K),  we  choose  these  methods  as  a  standard  for 


comparison.  The  k^  order  A-B  method  is  characterized  simply  by 

p(2)  =  zk_1  -  zk, 

so  it  is  the  member  of  P(k,  k)  where  all  the  parameters  designated 
as  free  in  section  1.3  have  been  set  to  zero.  Classically,  this 
has  been  regarded  as  a  desirable  property,  and  certainly  one  not 
possessed  by  the  members  of  P( k ,  k)  in  general.  It  was  thought 

4-U 

necessary,  for  a  general  k  order  method,  to  store  2k  values  between 
integration  steps  and  only  k+1  for  the  A-B  methods.  For  example  this 
notion  prevailed  in  limiting  the  search  of  Dill  and  Gear  [7]  for 
stiffly  stable  methods  to  this  type  of 'Inin  i  mum  storage  methods." 

j,  L, 

However,  Skeel  [38]  found  that  foe  P-C  methods  where  a  J,  order 
corrector  is  paired  with  a  specific  kth  order  predictor  it  is  necessary 
to  store  at  most  k+1  values,  regardless  of  zero  coefficients.  If  we 
apply  these  findings  to  codes  employing  predictors  alone  we  see  that 
any  k-step  method  needs  at  most  k  values  stored  between  integration 
steps.  This  removes  storage  economy  as  a  relevant  comparison. 

As  indicated  in  Table  4.2-1,  we  have  found  predictors  with 
stability  radii  roughly  three  times  that  of  the  A-B  method  of 
corresponding  order.  Another  valuable  perspective  is  gained  by 
noting  the  radius  of  any  order  A-B  methods  is  roughly  matched  by 
our  near-optimal  method  of  two  orders  higher.  Normally  an  increase 
in  order  means  more  accuracy  at  the  expense  of  computation  time 
induced  by  the  smaller  stepsize  lequired  of  reduced  higher  order 
stability  regions.  As  illustrated  later  in  this  section,  we  can  leave 
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the  stepsize  the  same  and  change  from  an  A-B  method  to  a  near- 
optimal  method  of  two  orders  higher. 

Table  4.2-1. 


Order 

Adams-Bash forth 
Radius  CjHl=Cp+1/a(TI 

kadi  us ' 

Near-Optimal 

Cp+1  Cp+l/a(1) 

4 

0.2146 

0.3486 

0.4706 

0.3856 

0.7711 

5 

0.1266 

0.3297 

0.3470 

0.3666 

1 .3609 

6 

0.0731 

0.3156 

0.2395 

0.3496 

7.4542 

7 

0.0412 

0.3042 

0.1595 

0.3259 

2.1067 

8 

0.0226 

0.2947 

0.0745 

0.3056 

1.4813 

9 

0.0121 

0.2870 

0.0405 

0.2957 

2.4933 

We  placed  no  control  on  the  expression  for  local  or  global 
truncation  error  constants.  However,  the  local  truncation  error 
constant  remained  approximately  the  seme  as  that  for  the  A-B  methods. 
The  constant  related  to  global  truncation  error,  Cp  i.-]/a(l ) ,  did 
increase  moderately.  This  is  a  Ju  ult  of  the  observation  that  as  we 
approach  optimally  stable  methods,  a  root  of  a(z)  goes  to  1.  It  also 
appears  there  is  no  appreciable  gain  in  stability  when  the  constant 
Cp+]/o(l)  is  increased  beyond  a  certain  point.  As  shown  by  example 
later  in  this  section,  the  increase  of  this  constant  for  these 
near-optimal  methods  is  not  so  severe  as  to  damage  their  utility. 
Another  study  addressing  the  question  of  how  far  these  stability 
regions  can  be  extended  for  fixed  error  constants  would  be  of  interest. 


Properties  of  these  near-optimal  predictors  are  given  in 
Table  4.2-1.  Rational  coefficients  for  these  methods  are  given  in 
Appendix  B.  Note  for  A-B  Methods,  c{l)  =  1  so  that  Cp+-|/o(l)  =  C  -j. 

To  illustrate  the  effectiveness  of  the  near-optimal  stability 
properties  given  in  Table  4.2-1  and  the  necessity  to  measure  relative 
stability  we  solve  two  example  problems  in  35-digit  precision  with 
A-B  methods  and  with  the  near-optimal  methods.  The  first  problem 
considered  is 

{/'  K  -16 y,  y{ 0)  =  1,  xe[0,  3]. 

He  use  the  A-B  5th  order  method  on  the  /'th  order  near-optimal  method, 
both  with  a  stepsize  of  h  =  0.01.  Thus  we  have  hX  =  0.16  which  is 
inside  both  methods'  region  of  absolute  stability.  However  it  is 
well  outside  the  A-B  5th  order  relative  stability  region  and  near 
the  boundary  of  relative  stability  for  the  7th  order  near-optimal 
method.  This  problem  gives  us  an  example  of  what  can  happen  when  we 
forget  about  relative  stability  in  the  left  half  plane.  If  we  had 
solved  this  problem  with  the  A-B  5th  order  method  and  desired  to 
follow  the  solution  to  this  problem  (possibly  a  component  of  a  larger 
system)  until  the  solution  fell  below  10-^  for  example,  we  would  find 
by  continuing  the  calculation  it  is  necessary  to  continue  to  step  825. 
Whereas  in  fact,  and  as  detected  by  the  relatively  stable  7th  order 
near-optimal  method,  we  could  terminate  calculation  after  step  150. 
in  this  case  and  others  like  it,  ignoring  the  relative  stability 
characteristics  of  a  method  can  be  a  very  costly  choice.  The 
relative  error  for  the  A-B  solution  in  Table  4.2-2  after  300  steps 


is  1 . 6000Qf*1 3  and  for  the  near-optimal  method  is  5.3610Q-03.  Note 
that  the  A-B  solution  is  tending  to  zero,  which  evidences  the 
methods'  absolute  stability  at  hX  =  -0.16,  but  not  nearly  so  fast 
as  the  solution  itself.  In  fact  it  tends  to  zero  so  slowly  that  the 
calculated  solution  values  are  essentially  meaningless. 

Table  4.2-2. 


Step 

No. 

True 

Solution 

A-B  5th  Order 
Calculated  Actual  Error 

7th  Order 
Calculated 

Near- Optimal 
Actual  Frror 

0 

0.  lOOOQi-Ol 

O.iOOOQ+Ol  0.0 

0. 1000Q+01 

0.0 

50 

0. 3355Q-03 

0. 3357Q-03  -0.1988Q-06 

0. 3354Q-03 

0.6280Q-07 

100 

0.11 25Q-06 

0. 3001Q-06  -0.1 875Q-06 

0.11 25Q-06 

0. 851 5Q--10 

150 

0. 3775Q-10 

0. 1108Q-06  -0. 1108Q-06 

0. 3769Q-10 

0.6116Q-1 3 

200 

0. 1266Q-1 3 

0. 6541Q-07  -0.6541Q-07 

0. 1263Q-1 3 

0.340 IQ- 16 

250 

0. 4248Q-1 7 

0. 3862Q-07  -0.3862Q-07 

0. 4231Q-1 7 

0. 1686Q-1 9 

300 

0.1425Q-20 

0.2280Q-07  -0.2280Q-07 

0. 1418Q-20 

0. 7640Q-23 

The  second  example  is  taken  to  show  simply  that  indeed  there 
are  values  of  hX  outside  the  A-B  methods’  stability  region  but  inside 
the  near-optimal  methods'  stability  region.  ’Besides  highlighting 
the  difference  in  the  stability  regions  this  example  provides  credence 
in  the  method  used  to  evaluate  stability.  The  problem  used  is 

</'  =  -y,  i/(0)  =  1 ,  xc[0,  10]. 

It  is  solved  by  the  7th  order  A-B  method  and  the  7th  order  near- 
optimal  method  given  in  Table  4.2-1.  A  stepsize  of  h  =  0.1  was  used 


iijjf  ,  * 


in  both  solutions  resulting  in  hX  -  -0.1  which,  according  to  the 
table,  is  outside  the  A-B  stability  region  but  inside  the  near 
optimal  stability  region.  The  results  are  given  below  in  Table  4.3-3. 

Table  4.2-3. 


Step 

No. 

True 

Sol  ution 

A-B  7th  Order 
Calculated  Actual  Error 

7  th  Order  [ 
Cal cul ated 

Joar-Ootimal 
Actual  Error 

0 

0.1 000Q+01 

0. 1000Q+01 

0.0 

0.1000Q-.-01 

0.0 

29 

0.8208Q-01 

0.8208Q-01 

0. 51 37Q-07 

0.8208Q-01 

0.4764Q-07 

50 

0.6738Q-02 

0. 6788Q-02 

-0.4977Q-04 

0.6738Q-02 

0.1 320Q-07 

75 

0 .5531Q-03 

-0.5459Q-01 

0. 5514Q-01 

0.5531Q-03 

0. 1S36Q-08 

100 

0.4540Q-04 

0.6109Q+02 

-0.6109Q-02 

0.4540Q-04 

0. 2304Q-09 

We  see  clearly  that  the  A-B  method  is  unstable  for  this  problem. 

Being  able  to  increase  the  order  by  two  without  reducing  the 
stepsize  more  than  offsets  the  larger  error  constants  possessed  by 
the  near-optimal  methods.  This  can  be  done  as  long  as  we  are  inside 
the  stability  region  and  the  stepsize  does  not  exceed  a  certain  size. 
For  example,  the  part  of  the  error  dependent  upon  method  and 
stepsize  for  the  7th  order  A-B  method  is  0.3042h^  and  for  the  9th 
order  near-optimal  method  is  2.4933h9.  Thus  for  h  <  0.3493  the  9th 
order  near-optimal  method  is  more  accurate  than  the  7th  order  A-B 
method.  We  compute  these  limiting  values  of  h  for  all  such 
comparisons  in  Table  4.2-4  below. 


Table  4.2-4 


Near-Optimal 

Li mi  ting 

Order 

h-Val ue 

We  note  all  these  limiting  values  of  h  are  well  outside 
normal  stepsize  ranges  and  that  for  h  less  than  these  values, 
increased  accuracy  will  be  realized  by  the  higher  order  near-optimal 
methods,  for  example,  in  comparing  the  6tn  order  A-B  method  with  the 
8th  order  near-optimal  method,  the  part  of  the  error  dependent  upon 
stepsize  and  error  constant  for  h  -  0.1  is  0 . 31 56Q-07  and  0.1481Q-08, 
respectively. 
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CHAPTER  V 


STIFFLY  STABLE  CORRECTOR  SEARCH 

5.1  Statement  of  and  Approach  to  the  Problem 

In  the  past  fifteen  years  much  research  has  been  devoted  to 
the  development  of  numerical  methods  for  obtaining  solutions  to 
problems  (1.1-1)  of  the  class  referred  to  do  stiff  problems.  These 
problems  and  the  difficulties  they  engender  were  described  in 
section  2.1.  Most  popular  codes  which  employ  LMM  for  the  solution  of 
stiff  problems  depend  exclusively  on  the  backwards  difference  formulae 
(BDr )  introduced  by  Gear.  These  methods  suffer  large  regions  of 
instability  in  the  left  half  plane  and  are  not  A^-stabie  for  orders 
greater  than  six.  This  suggests  much  improvement  is  possible. 

As  mentioned  in  Chapter  II,  Grigorieff  and  Scholl  [12]  and 
Kong  [25]  have  given  constructive  proofs  which  show  the  existence  of 
A(a)-stable  methods  of  arbitrarily  high  order  with  a  arbitrarily  close 
to  n/2.  However,  searches  for  methods  successful  in  the  solution  of 
stiff  problems  have  made  little  progress.  Kong  [25]  performed  a 
numerical  optimization  on  a  for  fixed  values  of  the  error  constant 
although  these  methods  have  little  practical  value  because  of  their 
extremely  small  regions  of  accuracy  about  the  origin.  In  a  series  of 
papers  [40,  15,  13],  Gupta  and  Wallace  found  methods  with  improved 
values  for  by  describing  LMM  in  terms  of  local  polynomial  approxi¬ 
mations.  The  methods  resulting  from  their  investigations  have 
larger  error  constants  which,  for  many  of  their  methods,  discount 
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the  practical  value  they  might  otherwise  have.  Dill  and  Gear  [7] 
tried  a  computer-aided  trial  and  error  approach  with  no  appreciable 
results.  None  of  these  investigations  have  resulted  in  methods  which 
perform  significantly  better  than  tiie  BDF. 

We  see  from  these  reports  there  are  several  properties  of  LMM 
which  appear  decisive  in  determining  the  performance  of  a  LMM  in  the 
solution  of  a  stiff  problem.  To  conduct  a  search  for  optimal  methods 
we  need  to  isolate  these  properties,  determine  a  useful  measure  for 
each  property,  and  define  our  object  function  in  terms  of  these 
measures.  The  region  of  absolute  stability  is  of  prime  interest 
because  cf  the  need  to  maintain  stability  for  the  components  resulting 
from  the  eigenvalues  with  negative  real  part  and  large  magnitude.  The 
angle  a  is  n  measure  of  the  absolute  stability  region  which  is  widely 
used.  It  is  a  good  indicator  of  the  number  of  problems  on  which  the 
LMM  may  be  used.  It  is  also  a  natural  measure  based  on  implementation 
of  LMM  in  that  stcpsize  is  the  factor  which  leads  to  difficulties  in 
the  solution  of  stiff  problems  with  LMM  possessing  finite  regions  of 
stability.  A  change  in  stepsize  moves  a  fixed  eigenvalue  radically 
toward  or  away  from  the  origin.  The  angle  a  indicates  whether  such 
radial  movement  may  intersect  regions  of  instability  for  a  given 
eigenvalue. 

Hie  region  of  accuracy  about  the  origin  is  important  because 
the  components  of  the  solution  resulting  from  eigenvalues  near  the 
origin  are  dominant  in  the  solution  and  their  accuracy  must  be  of 
concern.  A  natural  measure  for  this  region  is  the  radius  of  relative 
stability.  It  seems  clear  that  any  measure  for  this  region  must  admit 
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the  origin  as  an  interior  j;oint  unless  we  preclude  problems  with 
eigenvalues  of  positive  real  part.  Also,  as  vie  noted  with  the  problem 
solution  given  in  Table  4.2-2,  a  knowledge  of  the  relative  stability 
characteristics  is  essential  before  we  attach  meaning  to  the 
numerical  solution,  even  for  eigenvalues  with  negative  real  pa^t. 

Finally  the  error  term  is  a  decisive  property  since  as  vie  see  in  these 
previous  investigations,  it  can  become  prohibitively  large. 

These  considerations  led  to  the  definition  of  A(cx,  r)  stability 
given  in  section  2.3.  It  seems  of  interest  tner.  to  know  how  far  a 
and  r  can  be  extended  for  a  fixed  value  of  the  error  constant.  We 
report  the  results  of  such  an  investigation  in  this  chapter  and  find 
methods  nearly  optimal  with  regard  to  these  desired  properties. 

Further,  in  comparison  vie  find  these  near-optimal  methods  perform 
successfully  as  we  would  expect  from  the  stability  and  error  measure¬ 
ments  applied. 

Wc  choose  the  class  C(n,  n)  for  our  search  of  optimal  methods. 
Our  reasoning  for  this  choice  is  similar  to  that  u?ed  in  Chapter  IV 
when  we  chose  the  class  P(n,  n)  for  our  predictor  search.  Theorem 
4.1-1  limits  our  investigation  to  classes  C(n,  m)  where  n  £  m  +  1. 

Other  investigations,  such  as  those  previously  mentioned  were 
conducted  in  C(n,  n)  and  the  currently  popular  BDF  also  lie  in  this 
class.  These  methods  provide  good  comparisons  for  our  results.  In 
addition,  restricting  our  search  to  C(n,  n-1)  seems  to  be  an  unrewarding 
handicap.  We  did  make  several  optimization  runs  w'thin  the  class 
C(4,  5)  with  no  significant  improvement  over  results  we  obtained  from 
the  class  C(4,  4).  Upon  consideration  of  these  factors  we  conclude 
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C{n,  n)  is  the  better  class  for  our  search.  Further,  we  search 
for  A(a)-stnble  methods  and  since  then  ability  is  necessary,  by 
virtue  of  theorem  2.3-7  A^-stabili ty  is  also  necessary  provided  we 
allow  no  roots  of  a(z)  to  have  unit  magnitude.  Applying  the 
characterization  of  A  -stable  methods  from  theorem  2.3-6,  we  consider 
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only  the  subclass  of  C(n,  n)  for  which  all  roots  of  o(z)  are  of  less 
than  unit  magnitude. 

The  corrector  search  is  posed  as  an  optimization  problem 
similar  to  that  of  the  predictor  search  in  Chapter  IV.  The  program 
we  use  ib  basically  the  sai.s  with  five  notable  exceptions.  First  of 
all,  in  this  problem  we  have  two  desired  prooertios,  a  and  r,  and 
our  object  function  for  the  optimization  is  def’.ted  as  a  linear 
combination  of  these  two  properties.  The  second  difference  in  the 
two  programs  is  sizeable  and  derives  from  the  need  to  evaluate  g. 

The  remaining  changes  we  refer  to  are  those  necessary  for  fixing  the 
error  constant,  switching  to  correctors,  and  for  limiting  our  search 
to  A  -stable  methods. 

CO 

Defining  the  object  function  as  a  linear  combination  of  a  and 
r  increases  the  scope  of  the  optimizazion  problem  since  then  our 
interest  extends  to  the  effect  of  taking  different  linear  combinations. 
We  find  an  indicated  relationship  between  the  values  of  a  and  r 
corresponding  to  maximal  values  of  these  different  linear  combinations. 
The  two  properties  are  inversely  related  but  not  linearly  so.  For 
example,  the  following  relationships  are  indicated  to  exist  within 
C{4,  4). 
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radius  of  relative  stability 


Figure  5.1-1.  a  vs.  r. 


The  three  curves  result  from  fixing  the  error  constant. 


C,.+-j/o(l),  at  the  values  given  adjacent  to  the  curves.  Notice  the 


curves  seem  to  approach  the  limiting  a- value  of  Ii/2  for  larger  values 


of  the  radius  as  the  error  constraint  is  relaxed.  This  leads  one 


to  suspect  these  curves,  corresponding  to  values  of  /a ( 1 )  in 

(0,  «>),  fill  in  the  area  under  the  line  a  =  71/2  and  to  the  left  of 
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the  line  r  =  )'q,  where  is  the  largest  radius  of  any  method  in  the 
class  C(n,  n).  An  estimate  for  the  value  of  rQ  is  provided  by  the 
work  of  Thompson  and  Rodabaugh  [39],  where  they  report  candidates  for 
rQ  values  for  the  classes  C(n,  n-1),  4  <_  n  _<  9.  They  obtained  these 
values  by  a  numerical  optimization  with  no  constraint  on  the  error 
constant.  We  did  not  investigate  thoroughly  for  all  orders  and 
error  constants  the  portion  of  the  curve  defined  by  linear  coiibi nations 
in  which  the  radius  was  weighted  much  more  than  alpha.  The  methods 
in  this  portion  of  the  curve  are  expensive  to  find  and  of  little  value 
in  the  solution  of  stiff  problems  because  of  their  poor  absolute 
stability  properties.  Kong  [25]  considered  only  absolute  stability, 
and  therefore  his  results  define  the  upper  extent  of  the  curve  on  the 
alpha  axis.  The  optimization  runs  for  these  methods  are  inexpensive, 
however,  they  also  are  of  little  value  in  the  solution  of  stiff  or 
non-stiff  problems  because  of  their  poor  relative  stability  properties. 
For  orders  greater  than  four  we  concentrate  our  effort  in  the  region 
about  and  just  to  the  left  of  the  point  where  the  curve  starts 
dropping  off  most  rapidly.  This  region  of  the  curve  yields  methods 


which  pair  values  of  a  and  r  in  a  way  least  costly  to  either  property. 

The  second  change  is  that  necessary  for  evaluating  the  angle 
a  of  A(a)  stability.  We  use  the  stability  function 
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p(r)  we  !kv?  ,/x)  =  .,( z) 3  so  that  c  ( x ,'  =  ,-; { / '  and  j  not  cc.:c-  -:■; 

oursel \cs  with  the  pes  Mon  of  the  tr.’i  circle  be!...;  •.  ':-  rc..*  u>:': s . 

Cooke  [?  1  fir.  :s  ch.v  o  will  ho  pcfiti,'r  as  long  :'.:>  0(2}  ii  inject!  v„- 
.  I  analytic  in  i.iiG  c];?S‘  CT  t«2  C.;<  £!  io-“  Of  tht  I  mil  ci  re.':  2.  We 
guaraf. c.  ,:is  realm  of  i:u«lyf icily  !>>  ,\  -stability  aud  the  o"ly  ers's 
we  fine  in  M.  jecti  ,'i  ty  is  '  is  i"  creasing  '•R0  real  •'  x i s 

te  toe  lof.  of  the  v.  ‘  ■,  n.>J';i.  3  fro:-?  2  r.-thod  c  •'  no  \  in 

this  investig-Mm  since  ins.it  0  -  0.  Ii.se  boa-dary  cyiv-  r.:y  hfc.e 
po-itive  or  negative  orientation  (Are;  '.;(e_;  ’)  -  /.ig  ( q  { 1 ) ;  .  ?.\.) . 

Wo  incu.e  positive  criers L:-ti on  hy  ossi-grirv  tea  emulated  uaf.insry 
part  of  qfa)  t*>  be  positive.  For  evaluating  o  v:e  are  interested  or.iy 
in  those  boundary  points  with  negative  real  pari.  From  the  resulting 
points  of  interest  calculated*  we  find  the  one  w*th  the  largest  ratio. 


I  mag  (c.(e  n))/Se  (q(e  n)). 


Alpha  is  taken  as  the  arctangent  of  the  negative  of  this  largest 
ratio.  It  is  important  that  the  coefficients  cf  p  and  c  are 
associated  as  precisely  as  possible  to  a  member  of  C(n,  n).  We  found 
when  poor  precision  was  used  that  the  resulting  locus  calculated  was 
a  distortion  of  the  locus  we  intended  to  investigate.  This  was 
especially  true  for  v  near  the  origin  where  occasionally  the  loci 
approached  the ori gin  along  the  real  axis  instead  of  along  the 
imaginary  axis.  Within  C(n,  n)  for  example.  Theorem  2.2-4  guarantees 
this  cannot  happen.  This  distortion  resulted  in  meaningless  values 
for  a.  The  use  of  16-digit  precision  in  the  coefficients  representing 
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the  members  of  C(n,  n)  we  investigated  eliminated  this  problem.  The 
routine  described  here  to  evaluate  a  is  very  fast  and  its  cost  is  a 
relatively  minor  part  of  the  optimization. 


The  next  change  we  discuss  is  that  required  to  restrict  our 
investigation  to  the  subclass  of  C(n,  n)  corresponding  to  ?  fixed 


error  constant  C  .j/a(l)  =  Eg. 


The  relationship 


n-1 
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derived  from  requiring  =  0  in  the  system  (1.2-3),  coupled  with 
the  requirement  C  ^  =  Eg.  yields  an  expression  for  8 ,  explicitly  ir, 
terms  of  a. ,  ,  «n ^  and  EQ.  For  example  if  n  =  3.  we  get 

3.-,  =  C(9-24Eq)  -  (i  +  24S0)o1  -  4SEpa2]/?4.  (5.1-1) 


As  expected,  this  leaves  one  fewer  free  parameters  than  indicated  in 
the  discussion  of  section  1.3.  The  domain  of  our  object  function  for 
the  class  C(n,  n)  with  fixed  error  constant  is  then  Rn~J,  the  same 
as  for  the  class  P(n,  n)  with  a  floating  error  constant. 

The  last  two  changes  referred  to  earlier  were  those  required 
for  switching  to  correctors  and  insuring  A^-stability.  Working  with 
correctors  involves  inclusion  of  8  ^  in  calculations  such  as 
determining  the  coefficients  of  the  LMM  and  in  defining  the 
characteristic  polynomial.  A^-stability  is  determined  simply  by 
checking  the  magnitude  of  the  roots  of  o.  The  optimizer  is  constrained 
away  from  points  in  the  domain  with  roots  of  o  greater  than  one  by  use 
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of  the  penalty  function  described  in  [8].  The  use  of  this  constraint 
works  very  effectively  in  maintaining  A^-stabi lily. 

As  in  section  4.1,  to  guarantee  a  positive  radius  of 
relative  stability  and  hence  zero-stability  also,  we  consider  only 
those  methods  for  which  p{z)  =  (z-l)p^(z)  where  p^(z)  is  a  Schur 
polynomial.  Since  the  coefficients  of  a(z)  are  determined  by 
a] 5  ***  >  an_i»  and  £q>  -be  requirement  of  ^-stability  can  be  used 
to  further  refine  the  region  of  interest  for  our  investigation.  As 
an  example  we  examine  these  conditions  for  the  class  C(3,  3)  with 
fixed  error  constant  E^.  The  conditions  insuring  p<.  to  be  Schur  which 
follow  are  gained  f>om  an  application  of  Theorem  2.1-4. 

1  -  ot.j  >0 

1  -  a2 ,  >  0  (5.1-2) 

1  +  a.j  +  >  0 


Using  the  corollary  2.1-5  we  see  necessity  for  the  polynomial  o(z) 
to  be  Schur  is  similarly  provided  by  all  elements  of  the  product 
vector  (^(^2*  $y  &q>  &  -|)^  being  of  the  same  sign.  Using  the 
relationships  between  the  a.  and  3..  developed  in  section  1.3  and 
given  in  (5.1-1)  we  express  the  elements  of  the  producx 
M3 ($2 >  B-j »  3qj  3  1)^  in  terms  of  the  and  Eq.  Since  tv/o  of  these 
elemc.  <.s  are  the  same  as  two  of  the  expressions  in  (5.1-2)  required 
positive  for  p^  to  be  Schur,  all  elements  of  the  product  vector  must 
be  positive.  The  resulting  conditions  for  Aro-stabil  ity  are  given  in 
(5.1-3).  The  first  two  of  these  conditions  supplement  the  conditions 
for  zero-stability  given  in  (5.1-2). 
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-  (12EQ  +  1)  -  12EQa1  -  (24E0  -  1)  a,  >  0 

1  -  2  a-j  -  >  0 

(5.1-3) 

1  -  o2  >  0 

1  +  a-j  +  2  ag  >  0. 

The  first  inequality  of  (5.1-3)  is  the  only  express  ion  involving  Eg. 

This  class  of  methods  will  always  have  a  negative  error  constant  since 

when  Eq  >  0  the  first  and  third  inequalities  in  the  presence  of  the 

fourth  are  incompatible.  In  general,  since  3  ^  is  the  only  parameter 

introducing  terms  involving  Eg,  since  the  coefficients  of  3  ^  in 

each  of  the  expressions  for  3  ^ ,  g,  ,  6n  -j  respectively  define 

the  first  column  the  matrix  M  from  the  proof  of  theorem  2.1-4,  and 
?  n 

since  M*-2  In+^  we  will  have  the  error  constant  Eg  appearing  in  only 
one  of  the  inequalities  such  as  given  in  (5.1-3).  We  have  not  shown 
in  general  that  Eg  must  be  negative.  However,  as  in  (5.1-1),  equating 
the  expression  for  Cp+-j  to  Eg  •  o(l)  yields 

k-1 

6_-,  =  aQ  +  J  a.  a.  -  Ega(l)  (5.1-4) 

in  general  for  the  class  C(p,  k).  This  makes  it  clear  that  as  Eg  -» 

tiie  inequality  3_-j  approaches  the  inequality  o(l)  >  0,  which  is 

already  included  from  both  the  zero-stability  and  A^-stebil ity 

conditions.  Thus  as  we  would  expect,  for  large  negative  error  we 

gain  a  larger  region  of  interest.  This  leads  one  to  suspect,  and 

from  (5.1-2)  and  (5.1-3)  it  can  be  seen  for  C(3,  3),  that  if  a  given 

set  of  parameters  a,,  •••  ,  a  ,  determine  an  A  -stable  zero-stable 

i  n- 1 


method  for  any  error  constant  Eq  <  E^.  However  we  have  found  counter 
examples  to  this  conjecture  in  the  general  case. 

5 . 2  Results  and  Comparisons 

The  results  of  the  investigation  described  in  section  1  are 
given  here  for  orders  four  through  six.  The  properties  of  the  second 
order  trapezoidal  rule  and  the  third  order  BDF  suggest  little  or  no 
improvement  possible  below  fourth  order.  Since  the  BDF  are  members  of 
C(n,  n)  and  are  used  exclusively  in  popular  stiff  codes,  we  use  them 
for  comparison.  We  choose  a  selection  of  our  near-optimal  methods 
with  error  constants  C  ^/a(l)  of  -0.2  and  -0.5  for  fourth  order, 

-0.4  and  -0.8  for  fifth  order,  and  -0.9  for  sixth  order.  Their 
stability  properties  are  outlined  in  Table  5.2-1. 

Table  5.2-1.  Stability  Properties  of  Wear-Optimal  Stiff  Methods. 


Order 

Alpha 

Radius 

Wo(,) 

4 

1.481 

0.192 

-O.2C0 

1.414 

0.471 

1.377* 

0.650 

4 

1.535 

0.142 

-0.500 

1.511 

0.294 

1.445 

0.514 

5 

1.431* 

0.092 

-0.400 

1.338 

0.359 

1 . 259 

0.497 

5 

1.485 

0.077 

-0.800 

1.463 

0.155 

1.394 

0.463 

6 

1.321* 

0.121 

-0.900 

1.284 

0.299 

1.065 

0.435 

indicates  use  in  test  runs 

referred  to 

later  in  this 

section. 


For  comparison  we  give  correspond ing  information  for  the  BDF 
in  Table  5.2-2  below. 

Table  5.2-2.  Stability  Properties  of  the  BDF 


Order 

Alpha 

Radius 

C  . i /o( 1 ) 
pT  l 

4 

1.280 

0.484 

-0.200 

5 

0.905 

0.302 

-0.167 

6 

0.311 

0.130 

-0.143 

"fch 

The  kLn  order  BDF  has  an  error  constant  of  l/(k-H).  The 
severity  of  this  restriciion  on  the  error  constant  no  doubt  causes 
the  loss  of  Ag-stahility  for  the  higher  order  BDF.  Notice  for  each 
order,  the  near-optimal  methods  have  considerably  greater  regions  of 
relative  and  absolute  stability.  The  effect  of  the  larger  error 
constants  is  nullified  by  the  capability  to  increase  the  order  with  no 
accompanying  restriction  on  alpha.  For  example,  if  we  compare  the 
near-optimal  sixth  order  A(1.284,  0. 299)-stable  method  from  Table  5.2-1 
with  the  fourth  order  BDF,  it  is  the  case  that  the  problem  independent 
part  of  the  error  term  for  the  near-optimal  method  will  be  less  than 
that  of  the  BDF  for  stepsizes  less  than  0.471.  This  particular 
compar’son  does  not  involve  a  sacrifice  in  the  permissible  values  of 
alpha,  and  in  other  similar  comparisons  we  could  realize  again  in  the 
permissible  values  of  alpha. 

We  see  then  an  increase  in  the  error  constant  accompanying 
methods  with  higher  values  of  alpha  is  a  negative  aspect  that  can  be 
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controlled.  However,  an  inferior  value  of  alpha,  such  as  possessed  by 
the  fifth  and  sixth  order  BDF,  severely  disables  a  method.  We 
illustrate  this  fact  by  solving  two  systems  of  the  form  tj'{l)  -  i\y{ t) 
where  A  is  a  constant  2x2  complex  matrix  and  te[0,  5].  In  both 
problems  the  initial  values  are  determined  by  setting  both  constants 
in  the  general  solution  to  1.0.  Denote  the  problems  PI  and  P2  and  let 
them  be  defined  respectively  by  the  matrices  A1  and  A2  where 


A1 


-1  100 
0  -100  +  373i 


A2 


-1  100 
0  -100  +  250i 


(5.2-1) 


We  use  the  constant  stepsize  of  h  -  0.005  throughout  the  interval 
[0,  5].  The  calculated  values  are  given  in  Appendix  F.  The  near- 
optimal  methods  used  are  identified  by  an  asterisk  in  Table  5.2-1. 

We  summarize  the  results  of  these  computations  in  Table  5.2-3  below. 
All  calculations  were  done  in  35 -digit  precision. 

Table  5.2-3.  Summary  of  Tests  on  Problems  PI  and  P2. 


Order  of 


Method 

Method  Used 

Problem 

Res ul ts 

4 

A(1 .377,  0.650)-Slable 

PI 

Stable 

*  ■. 

t  l 

i  ^ 

4 

BCF 

PI 

Unstable 

i 

5 

A( 1 .431 ,  C.092)-Stable 

P2 

Stable 

5 

BDF 

P2 

Unstable 

\  * 

6 

A(1 .321 ,  0.121) -Stable 

P2 

Stable 

|  l 

6 

BDF 

P2 

Unstable 

\  f' 

Both  PI  and  B2  are  complicated  by  the  problem  of  stiffness. 
Notice  the  near-optimal  sixth  order  method  solved  the  problem  P2 
accurately  whereas  the  fifth  order  BDF  was  unstable. 

The  difference  in  the  relative  stability  radius  is  important 
as  we  see  in  the  following  example.  Consider  the  problem 

tj'  =  -8 tj,  </(0)  =  1,  xe[0,  10].  (5.2-2) 

We  solve  this  problem  with  the  fourth  order  A( 1 . 377 „  0.65G)-stable 
method  and  the  fourth  order  BDF  using  a  constant  stepsize  of  h  =  0.1. 
The  results  are  given  in  Table  5.2-4  which  follows. 


Table  5.2-4.  Solution  of  Problem  (5.2-2). 


Step 

No. 

True 

Sol  ution 

Fourth  Order  BDF 

A( 1 . 377,  C. 

.650) -Stable 

Calculated 

Actual  Err>r 

Calculated 

Actual  Error 

0 

0. 1000Q+01 

0. 1000Q+01 

0.0 

0. 1000Q+01 

0.0 

25 

0.205 IQ-03 

0. 2070Q-05 

-0.2050Q-OS 

-0.2179Q-03 

0.424CQ-08 

50 

0.4248Q-17 

0. 3318Q-12 

-0. 3318Q-12 

0 . 4 1 31 Q- 1 7 

0.117SQ-18 

The  previous  test  problems  PI  and  P2  demonstrated  the  need  for 
methods  with  higher  values  of  alpha  whereas  this  problem  demonstrates 
the  possible  effect  of  smaller  regions  of  accuracy  about  the  origin. 
The  relative  error  after  50  steps  in  0.02773  and  0.7811.10^  for  the 
near-optimal  and  BDF  methods,  respectively. 

We  close  with  one  more  demonstration  which  accents  the  value 
of  these  near-optimal  methods  as  multi-purpose  methods.  Lambert  [28, 


n  m  *L  ^  '  -M-VA.Vikn.V1 


ii.  ine  is  ■  of  stiff  rethrc!,  fr-r 


p.  ^ 79  j  indie  a ccs 

C  f  1  ul  'v*  1  V-S  t 

c  i  ;>  i ».  ine 

is  ■  of  stif 

f  re  tin  re! ,  if.- 

non-  s  .1  rf  problems .  Vi'  1.  • 

i)  la;  go  .  eepe 

c.  vf  accur 

.  .■  a  to  ‘t  ri.e 

origin  tnese  ne.-r 

-(.psivu  !\  / 

.a.  r)  -  stable  :• 

etKds  sr.ou 

i o  be  ideal  1^ 

suited  icr  this  p 

urposo.  f.^n 

■  idol  the  niv-o. 

L ' 

i  1  =  -.y 

'!{  ')  ~  1  x-  "C 

-  i]. 

(5.2 

We  solve  this  problem  with  th 

e  four  c'n  order 

AcLm  vu.i  netfiocl  and 

the  fifth  order  A{1.259,  j. 437)- stable  msiho 

■d  from  Tabl 

e  5.2  1.  ihe 

results  given  in 

Tabic  5.2-5  ■ 

sere  obt-  b.ed  us 'by,  a  roes 

tant  stepsize 

of  0.01. 

Tabic  5.2-5.  Solution  of  P r 

1  1  .  ! r.  p.  ->\ 

o  i  v.  I.  (  c  »  C  j 

Step  True 

No.  Solution 

4t'n  prde 
Ca  I  cti  i  ated 

Act  toil  t.rrcr 

5  th  Or^or 
Cal  cola  ,od 

!\ecr;diytiy’a  \ 
nC t.rai  hr* or 

0  0.100QQ+01 

0. ICOOQ+Ol 

0.0 

0.1000'/-01 

0.0 

25  0.7788 

0.7788 

0.45570-10 

0.7788 

0. 52 71 Q- 11 

50  0.6065 

0. 6065 

0. 7532b- 10 

0.6065 

0.1 114Q-1C 

75  0.4724 

0.4724 

0.9046^-10 

0.4724 

0. 1355Q-10 

100  C . 3679 

0.367s 

0 .9491 Q- 10 

0.3679 

0  1436Q-10 

As  expo- ted,  we  see  the  large  difference  in  error  constants 
for  these  two  methods  (-0.02639  and  -0.4  for  the  Adams-Mcul ton  and 
near-optimal  mc-.hods,  respectively)  is  more  than  offset  fo»*  t.ms 
stepsise  by  im;ng  a  method  of  one  order  higher.  An  alternative 
method  of  comps  sating  or  the  larger  error  constant  is  by  decreasing 
the  stepsize  bu  retaining  the  same  order.  We  do  this  with  the 
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fourth  order  A(1.377,  0. 650)-stable  method  and  a  stepsize  of  0.005. 
The  results  are  given  in  Table  5.2-6. 


Table  5.2-6. 

Solution  of  (5.2-3) 

by  Fourth  Order 

Nea r-Opti mal  Method . 

Step 

True  Solution 

Cal culated 

Actual  Error 

0 

0 . 1 00DQ+01 

0. 1000Q+01 

0.0 

50 

0.7780 

0.7783 

0.2269  Q-10 

100 

0.6065 

0 . 6065 

0.3580  Q-10 

150 

0.4724 

0.4724 

0.4356  0-10 

200 

0.3679 

0.3679 

0.4553  Q-10 

We  see  the  error  for  this  fourth  order  solution  is 
approximately  one-lvalf  that  of  the  fourth  order  Adams -Moulton  solution 
given  in  Table  5.2-5.  i-rom  an  examination  of  the  problem  independent 
part  of  the  methods'  error  terms,  we  would  expect  more  accuracy  from 
the  near-optimal  method  as  long  as  the  ratio  of  the  stepsizes  was  less 
than  0.6. 

In  this  investigation  we  have,  found  a  large  number  of  near- 
optimal  methods  and  it  is  neither  piactical  nor  necessary  to  list  the 
coefficients  of  all  these  methods.  For  those  methods  listed  in 
Table  5.2-1,  we  give  the  rational  coefficients  in  Appendix  C.  In 
Appendix  D  we  chart  the  properties  of  many  selected  methods  which 
indicate  the  relationship  between  alpha  and  the  radius  discussed  in 
the  previous  section.  A  relationship  similar  to  that  is  indicated  to 
exists  between  alpha  and  the  radius  when  instead  of  fixing  the  error 


For  the 


constant  C^+^/o(1),  we  fix  the  error  constant  C^/a  -| 
class  C(4,  4)  we  chart  the  latter  relationships  in  Appendix  t.  It 
is  interesting  to  note  how  well-defined  these  curves  are,  evidently 
indicative  of  the  geometry  of  the  range  for  the  function  with  which 
we  work. 

There  remains  questions  of  interest  in  this  area  which  have  not 
been  investigated.  For  example,  as  pertains  to  stiff  systems,  we 
clearly  need  to  determine  whether  results  as  rewarding  as  we  have  found 
here  await  our  investigation  into  higher  orders.  Ir  addition,  there 
may  exist  unconventional  methods  which  could  be  successfully  applied 
to  large  classes  of  stiff  problems.  Tiieie  currently  exists  many 
computational  difficulties  which  must  be  dealt  with  in  stiff  codes 
employing  LMM.  For  example,  since  the  implicit  relation  is  solved  by 
a  form  of  hew ton  iteration,  no  small  difficulty  is  presented  by  the 
need  to  calculate  me  Jacobian,  especially  for  large  systems.  Even  so 
we  find  LMM  are  the  most  popular  methods  used  to  solve  stiff  systems. 
This  indicates  great  potential  exists  for  a  different  approach.  As 
pertains  to  non-stiff  methods,  there  is  an  open  question  ol  how  best  to 
pair  predictors  and  correctors  in  P-C  methods.  Do  near-optimal 
predictors  and  near-optimal  corrector  pair  to  yield  near-optimal 
P-C  algorithms?  Also  as  mentioned  in  Chapter  IV,  the  question  of  how 
far  the  stability  regions  of  predictors  can  be  extended  for  fixed  error 
constants  remains  open. 
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In  the  matrices  given  below  we  list  only  the  numerator  of 

the  entries.  All  entries  of  the  matrix  for  C(n,  n)  have  the  common 

denominator  D  . 

n 


"  55  9  8  9  -So" 

-59  1  9  32  27  144 

37  -5  8  27  -96 


-9  1  0  9  24 


C(5,  5),  Dg  =  720 


[1901 

251 

232 

243 

224 

-360 (F 

-2774 

646 

992 

918 

1024 

7200 

2616 

-264 

192 

648 

384 

-7200 

-1274 

106 

32 

378 

1024 

3600 

_  251 

-19 

-8 

-27 

224 

-720 

79 


4277 

475 

448 

459 

CO 

475 

-8640" 

-7923 

1427 

2064 

1971 

2048 

1875 

21600 

9982 

-798 

224 

1026 

768 

1250 

-28900 

-7298 

482 

224 

1025 

2048 

1250 

21600 

2877 

-173 

-96 

-189 

448 

1875 

8640 

1 

cr> 

27 

16 

27 

0 

475 

1440_ 

C(7,  7) ,  P7  =  60480 


198721 

19087 

18224 

18495 

18304 

L.O 

o 

co 

17712 

-423360"! 

-447288 

65112 

90240 

87480 

89083 

87000 

93312 

1270080  j 

705549 

-46461 

528 

31347 

24576 

31875 

11654 

-2116800 

-688256 

37504 

21248 

58752 

96256 

80000 

1 1 7504 

2116800 

407139 

-20211 

-12912 

-19683 

11136 

58125 

11664 

-1270080 

-134472 

6312 

4224 

5832 

3072 

28200 

93312 

423360 

19087  -863  -592 


-783  -512  -1375  17712 


-'"9489 


80 


C(8,  8), 


D9  120960 


"  434241 

36799 

-1152169 

1 39849 

2183877 

-121797 

-2664477 

123133 

2102243 

-88547 

-"*041723 

41499 

295767 

-11351 

-36799 

1375 

35775 

35424 

183625 

186624 

34875 

23328 

2081 25 

235008 

68125 

23328 

85275 

186624 

-12735 

35424 

1375 

0 

35424 

35775 

35584 

187648 

183465 

185344 

-20448 

37179 

27648 

78336 

160029 

228352 

-61664 

-81891 

-13568 

29952 

37179 

27648 

8352 

-10071 

-8192 

1024 

1215 

1024 

36799 

-967660 

175273 

3386880 

64827 

-6773760 

146461 

8457200 

146461 

-6773760 

64827 

3386880 

175273 

-967680 

36799 

120950 

l\ 


C(9 ,  9),  Do  =  3628300 


14  097  247 
-43  125  206 
95  476  786 


1  070  017 
4  467  094 
-4  604  594 


1  036  064 
5  842  688 
•1  359  808 


1  043  361 
5  743  062 
278  478 


1  040  128 
5  779  456 
62  464 


139 

855 

262 

5 

595 

358 

3 

842 

816 

6 

474 

654 

8 

384 

512 

137 

968 

480 

-5 

033 

120 

-3 

715 

840 

-4 

548 

960 

-2 

324 

480 

-91 

172 

642 

3 

146 

338 

2 

391 

296 

2 

789 

154 

2 

363 

392 

38 

833 

486 

-1 

291 

214 

- 

996 

928 

-1 

139 

022 

-1 

012 

736 

-9 

664 

105 

312 

874 

243 

968 

275 

562 

249 

856 

1 

070 

017 

-33 

953 

-26 

656 

-29 

889 

-27 

392 

1 

042 

625 

1 

039 

392 

1 

046 

689 

1 

012 

736 

-32 

659 

200 

5 

753 

750 

5 

785 

344 

5 

716 

438 

6 

029 

312 

130 

636 

800 

188 

7b0 

46 

656 

340 

942 

- 

950 

272 

-340 

819 

200 

7 

958 

750 

8 

356 

608 

7 

601 

566 

10 

747 

904 

457 

228 

800 

- 

100 

000 

- 

933 

120 

384 

160 

-4 

648 

960 

-457 

228 

800 

4 

273 

250 

6 

905 

088 

5 

152 

546 

10 

747 

904 

304 

819 

200 

-1 

228 

750 

409 

536 

3 

654 

322 

- 

950 

272 

-130 

636 

800 

286 

2r0 

186 

624 

1 

562 

218 

6 

029 

312 

32 

659 

200 

-30 

625 

-23 

328 

-57 

281 

1 

012 

736 

-3 

628 

800 

nflfiHfi 


The  following  coefficients  are  given  for  the  methods 
referred  to  in  Table  4.2-1.  The  numerators  of  the  coefficients  are 
listed  in  the  tables  below  and  each  order  has  a  common  denominator 
denoted  by  D. 

Table  B-l. 


Coeffi  cient 

4th  Order 

5  th  Order 

6  th  Order 

*. 

a-l 

-300 

-720  000  000 

-14  400  000 

3 

■r 

an 

LSI 

1  705  707  360 

40  006  080 

< 

0 

i 

ai 

-462 

-1  697  984  640 

-39  241  440 

A 

i 

i 

a2 

201 

987  585  120 

17  628  480 

% 

f 

i 

i 

a3 

-30 

-297  991  440 

-8  373  600 

t 

/* 

a4 

22  683  600 

6  521  760 

~ 

a5 

-2  141  280 

- 

eo 

570 

1  533  770  509 

33  963  773 

s 

el 

-869 

-3  284  474  686 

-97  823  787 

b2 

592 

3  245  856  024 

119  961  838 

63 

-143 

-1  604  272  786 

-81  922  322 

84 

303  066  559 

32  649  093 

65 

-6  153  235 

t 

D 

300 

720  000  000 

14  400  000 

t 

wmj~rr 


i  m 


Coefficient 


7th  Order 


8th  Order 


9  th  Order 


- 

■604 

CO 

o 

o 

GOO 

-1 

209 

600 

000 

-36 

288 

000 

000 

1 

479 

522 

240 

2 

352 

430 

080 

68 

929 

056 

000 

1 

562 

984 

640 

-2 

067 

206 

400 

-57 

556 

396 

800 

1 

112 

045 

76  D 

1 

67/ 

473 

280 

44 

478 

201 

600 

- 

•599 

477 

760 

- 

•895 

950 

720 

-29 

756 

160 

000 

252 

383 

040 

456 

624 

000 

25 

358 

054 

400 

-76 

688 

640 

• 

•411 

022 

080 

-13 

684 

204 

000 

0 

12 

096 

000 

-2 

474 

841 

600 

85 

155 

840 

232 

243 

200 

762 

O/’O 

000 

1 

698 

536 

399 

3 

982 

004 

743 

131 

051 

87y 

512 

5 

101 

980 

072 

-12 

450 

561 

759 

-461 

208 

069 

296 

8 

017 

316 

643 

23 

395 

159 

971 

1  008 

365 

805 

136 

7 

743 

202 

432 

-28 

566 

654 

731 

-1  467 

668 

483 

312 

4 

524 

155 

853 

22 

109 

841 

909 

1  435 

089 

863 

680 

1 

510 

919 

256 

-11 

107 

55i 

069 

-956 

912 

924 

432 

209 

655 

425 

3 

248 

488 

993 

403 

247 

072 

176 

-361 

187 

577 

-98 

967 

442 

256 

n 

306 

055 

592 

604 

800 

000 

1 

209 

600 

000 

36 

288 

000 

000 

<fe~2~&3«.3sSM 


*-  -  -  '  . 


The  following  coefficients  are  given  for  the  methods  referred 
to  in  Table  5.2-1.  The  numerators  of  the  coefficients  are  listed  in 


the  table  below.  The  coefficients  of  each  method  have  a  common 
denominator  denoted  by  D. 


Table  C-l. 

Fourth  Order,  C^/o(l) 

=  -0.200. 

A(1.481,  0192)- 
Coefficients  '  Stable 

A(1.'H4,  0.471)- 
Stable 

A(1 .377,  0.650)- 
St.  fie 

a-l 

-30  000  000 

-900 

-240  000 

°0 

60  267  000 

1800 

464  400 

al 

-38  319  000 

-1260 

-306  600 

°2 

1  743  O'O 

360 

87  283 

°3 

6  303  000 

0 

-5  088 

P-1 

14  032  375 

415 

10''  512 

eo 

1  137  750 

50 

24  165 

*1 

-10  470  000 

-240 

-64  993 

5  782  250 

no 

19  199 

^3 

3  551  625 

25 

4  829 

D 

30  000  000 

900 

240  000 

37 


Table  C-2.  Fourth  Order,  C./o{1)  =  -0.500. 

o 


A(1  .535,  0.142}-  A(  1 . 511  ,  0.294)-  A(1.445,  0.5U)- 

Coefficients  Stable  Stable  Stable 


“-1 

-900 

000 

000 

-900 

000 

-90 

000 

000 

“o 

2  312 

370 

000 

2  448 

009 

202 

104 

000 

-2  057 

940 

000 

-2  520 

000 

-160 

425 

COO 

“2 

630 

900 

000 

1  170 

000 

53 

032 

000 

C3 

14 

670 

G0C 

-193 

000 

-5 

661 

000 

*-1 

434 

431 

625 

437 

675 

45 

458 

925 

6o 

-231 

1  I'd 

750 

-31 7 

45C 

-20 

035 

950 

61 

-377 

409 

000 

-244 

200 

-9 

652 

200 

B2 

305 

311 

7  50 

329 

050 

1 

801 

550 

63 

16 

635 

375 

-79 

075 

2 

941 

675 

D 

900 

000 

000 

900 

000 

90 

000 

OCO 

88 


Table  C-3  Fifth  Order,  Cg/o(l)  =  -0.400. 


A { 1 .431 ,  0.092)-  A ( 1 .333,  0.399)-  A(1.259,  0.514)- 

Coefficienti  Stable  Stable  Stable 


a-l 

-9 

000 

000 

-9 

oeo 

000 

-720 

000 

ao 

27 

720 

000 

24 

120 

000 

1  800 

000 

al 

-35 

100 

COO 

-26 

100 

000 

-1  728 

000 

a2 

22 

500 

000 

14 

403 

000 

792 

000 

a3 

-7 

200 

000 

-3 

600 

O^o 

-144 

000 

a4 

1 

080 

000 

180 

000 

0 

B-1 

4 

159 

075 

4 

013 

625 

321 

408 

Bo 

-4 

117 

125 

-1 

923 

375 

-101 

£•40 

B1 

-1 

212 

750 

-2 

450 

250 

-202 

720 

e2 

4 

067 

250 

2 

829 

750 

117 

120 

B3 

-1 

537 

125 

• 

-693 

375 

38 

240 

B4 

259 

875 

- 

-156 

375 

-28 

208 

D 

9 

000 

000 

9 

000 

000 

720 

000 

E 


f 

| 

f 


Table  C-4.  fifth  Order,  C-/o{l) 


-0.800. 


Coe  ffi  dents 

A(1 . 

485,  0.077)- 
S table 

A{1 .46?,  0.155)- 
StJile 

A( 1.394,  0.463)- 
Stable 

a-l 

-9 

000 

000 

-14  400 

-9 

000  000 

a0 

28 

800 

coo 

44  640 

27 

352  800 

°1 

-36 

QfiO 

000 

-54  720 

-33 

753  600 

°2 

20 

700 

000 

31  680 

21 

635  900 

a3 

-4 

Vf 

* 

000 

-7  200 

-7 

236  000 

a4 

0 

0 

9S9  500 

8-i 

4 

210 

125 

6  711 

4 

1 '5  335 

p’0 

-4 

726 

875 

-6  833 

-5 

895  235 

6, 

t 

-2 

on 

250 

-2  998 

-1 

489  no 

62 

5 

iso 

750 

7  242 

3 

4t  690 

83 

-1 

576 

875 

-2  353 

-1 

577  785 

B4 

-191 

875 

-329 

223  905 

D 

S 

COO 

000 

14  400 

9 

000  000 

Table  C-5.  Sixth  Order,  C.;/c(l)  =  -0.900. 


Coefficient 


A(  1 . 321 ,  0.121)- 
Stable 


A(1.284,  0.299)-  A(1.065,  0.435)- 

S table  Stable 


ft-l 

-14 

400 

000 

-604 

800 

000 

- 

-604 

800 

ao 

49 

824 

000 

2 

04  7 

065 

SCO 

1 

959 

552 

al 

-72 

000 

000 

-2 

903 

040 

003 

-2 

721 

500 

i 

C2 

53 

280 

000 

2 

177 

280 

OCO 

2 

115 

800 

\ 

“3 

-18 

720 

000 

-846 

720 

000 

-967 

680 

A 

a4 

1 

440 

000 

117 

996 

489 

241 

920 

5» 

f 

T 

°5 

57C 

000 

12 

216 

SfO 

-24 

192 

% 

P-1 

6 

553 

520 

254 

974 

986 

274 

447 

3o 

-9 

114 

120 

-302 

886 

COO 

-330 

69  j 

* 

*1 

1 

415 

800 

-34 

322 

494 

114 

147 

; 

h 

4 

667 

600 

310 

644 

72  5 

-40 

448 

*3 

-1 

129 

200 

-148 

424 

994 

103 

437 

-1 

858 

120 

-14 

144 

400 

-81 

720 

e5 

894 

520 

23 

388 

886 

21 

313 

i 

D 

14 

400 

000 

604 

800 

000 

604 

800 

; 

APPENDIX  E 


RELATIONSHIP  OF  ALPHA  VS.  RADIUS  FOR  FIXED  LOCAL 
TRUNCATION  ERROR  CONSTANTS 


Figure  E-1.  Relationship  for  C(4,  4). 


APPENDIX  F 


SOLUTION  OF  TEST  PROBLEMS  PI  AND  P2 


Calculations  below  were  carried  out  on  an  Amdahl  370/V7  in 
35-digit  precision.  Numbers  less  than  approximately  10-^  in  magnitude 
are  represented  as  0.0. 


Table  F-l.  True  Solution  of  Problem  PI. 


Step 

(h=0.005) 

First  Component 

Second  Conronent 

Real 

Imaginary 

Real 

Imaginary 

0 

0.2000Q+01 

0.0 

-0.9900Q+00 

0.3730Q*01 

200 

0. 3679Q+00 

0.2794Q-43 

-0.7983Q-43 

-0.11S3Q-42 

400 

0.1353Q+00 

0.0 

0.0 

0.0 

600 

0.4979Q-01 

0.0 

0.0 

0.0 

800 

0.1832Q-01 

0.0 

0.0 

0.0 

1000 

0.6738Q-02 

0.0 

0.0 

0.0 

Table  F-2. 

Calculated  Solution  of  PI  by  Fourth  Order  A(l. 
Stable  Method. 

377,  0.650)-- 

Step 

Fir'S t  Component 

Second  Component 

(h=0.005) 

Real 

Imaginary 

Real 

Imaginary 

0 

0. 20000*01 

0.0 

-0.9900Q+00 

0.3730Q+01 

200 

0.3679Q+00 

-0.1916Q-05 

0.6626Q-05 

G.3862Q-05 

400 

0.1353Q+00 

0.4772Q-11 

-0.1198Q-10 

-0.2664Q-10 

600 

0.4979Q-01 

-0 . 1 773Q-*.  7 

-0. 2188Q-1 6 

0.1091 Q-15 

800 

0.1832Q-01 

-0.5826Q-22 

0. 3095Q-21 

-0.2697Q-21 

1000 

0.6738Q-02 

0.3869Q-27 

-0. 1601 Q- 26 

0.2124Q-27 

1000 


Table  F-3.  Actual  Global  Error  of  Calculated  Solution  of  PI  by 
Fourth  Order  A(1.377,  0.650)-Stabla  Method. 


Step 

First 

Component 

Second  Component 

Real 

Imaginary 

Real 

Imaginary 

0 

-0.5267Q-06 

0.1916Q-05 

-0.6626Q-05 

-0.3862Q-05 

0.3970Q-10 

-0.4772Q-11 

0-1 19SQ-10 

0.2664Q— 10 

600 

0.1872Q-10 

0. 1773Q-17 

0.2118Q-16 

-0.1019Q-15 

0.9199Q-11 

0.5826Q-22 

-0.3095Q-21 

0.2897Q-21 

0.4234Q-1 1 

-0. 3869Q-27 

0.1601Q-26 

-0.2124Q-27 

i 

i 

l 

l 

l 


Table  F-4. 

Calculated  Solution  of  PI  by  Fourth  Order  BDF. 

Step 

First  Component 

Second  Component 

(h=0.005) 

Real  Imaginary 

Real 

Imaginary 

0 

0.2000Q+01  0.0 

-0.9900Q+00 

0.3730Q+01 

200 

0.9475Q+01  0.596aq+01 

-0.31 28Q+02 

0.28060+02 

400 

0.2122Q+03  -0.13480+03 

0.2928Q+03 

0.9245Q+03 

600 

-0.8181Q+03  -0.5742Q+04 

0.2223Q+05 

-0.2633Q+04 

800 

-0.1285Q+06  -0. 3738Q+05 

0.2669Q-+06 

-C-.4424Q-+06 

1000 

-0.20080+07  0. 2348Q+07 

-0.6769Q+07 

-0.9815Q+07 

1000 


A..'?*? 


Table  F-5.  Actual  Global  Error  of  Calculated  Solution  of  PI  by  Fourth 
Order  BDF. 


Table  F-6.  True  Solution  of  Problem  P2. 


Step 

h=0.005) 


First  Component 
Real  Imaaina 


Second  Component 
al  Iinaginar 


TiTil 


0.200Q+01 

0.0 

-0.9900Q+00 

0.2500Q+01 

0.3679Q+00 

—0. 3610Q-43 

0.8139Q-43 

0.5816Q-43 

0.1353Q+00 

0.0 

0.0 

0.0 

0.4979Q-01 

0.0 

0.0 

0.0 

0.1832Q-01 

0.0 

0.0 

0.0 

0. 67380-02 

0.0 

0.0 

0.0 
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Table  F-7. 

Calculated  Solution  of  P2  by  F 
Stable  Method. 

ifth  Order  A{  1.431 

,  0.092)- 

Step 

First  Component 

Second  Component 

(h=0.005) 

Real  Imaginary 

Real 

Imaginary 

0 

0.2000Q+P1  0.0 

-0.9900Q+00 

0.25C0Q+01 

200 

0. 3679Q+00  -0.5067Q-12 

-0.8977Q-12 

0.5967Q-11 

400 

0.1353Q+00  -0.4956Q-23 

0.7781Q-22 

-0.1604Q-21 

600 

0.4979Q-01  -0.17f.7Q- 31 

-0.4509Q-32 

0.4382Q-32 

800 

G. 18320-01  -0. 6335 Q- 32 

0.2243Q-42 

-0.1293Q-42 

1000 

0.6738Q-02  -0.1333Q-32 

-0.9632Q-53 

0.21 67Q-53 

Table  F-8. 

Actual  Global  Error  of  Calculated  Solution  of  P2 

!  by  Fifth 

Order  A(1.431,  0.092) -Stable  Method. 

First  Component 

Second  Component 

Step 

Real  Imaginary 

Real 

Imaginary 

0 

0.0  0.0 

0.0 

0.0 

200 

-0.2413Q-11  0.5067Q-12 

0.8977Q-12 

-0.5967Q-11 

400 

-0.1687Q-12  0.4946Q-23 

-0.7781Q-22 

0.1 604Q-21 

600 

-0.9338Q-13  0. 1757Q-31 

0.4609Q-32 

-0.4982Q-32 

800  -0.4588Q-13  0.G335Q-32  -0.2243Q-42  0.1293Q-42 

-0.2112Q-13  0.1333Q-32  0.9632Q-53  -0.2167Q-53 


1  I 

I 

5  „ 


r  '-m 

-3^ 

-v**( 

r  -M 


■*3S 

-5KI 


4ttl 


1000 


a 


Table  F-9.  Calculated  Solution  of  P2  by  Fifth  Order  BDF. 


Step 

(h=0.G05) 

First  Component 

Second  Component 

Real 

Imaginary 

Real 

Imaginary 

0 

0.2009Q+01 

0.0 

-0.9900Q+00 

0.2500Q+01 

0.93050+03 

0.2558Q+04 

-0.7321Q+04 

-0.1914Q+03 

0.1696Q+08 

0.6429Q+C8 

-0.1775Q+09 

-0.2124Q+08 

0.2666Q+12 

0.1601Q+13 

-0.4267Q+13 

-0.9187Q+12 

800 

0.2852Q+16 

0.39530+17 

-0  - 1 01 6Q-+1 8 

-0.3200Q+17 

-0. 2026Q+-20 

0.9673Q+21 

-0.2398Q+22 

-0.1008Q+22 

5 

i 

i 

j 


Table-  F-10.  Actual  Global  Error  of  Calculated  Solution  of 
Fifth  Order  BDF. 

■  P2  by 

First  Component 

Second 

Component 

Step 

Real  Imaginary 

Real 

Imaginary 

0 

0.0  0.0 

0.0 

0.0 

200 

-0.9362Q+03  -0.2558Q+04 

0.7321 Q+04 

0.1914Q+03 

400 

-0.1696Q+08  -0.6429Q+08 

*  0.1775Q+09 

0.2124Q+08 

600 

-0.2666Q+12  -0.1601Q+13 

0.4267Q+13 

0.9187Q+12 

800 

-0.2852Q+15  -0. 3953Q+17 

0. 1016Q+18 

0.3200Q+17 

1000 

0.2026Q+20  -0.9673Q+21 

0.2398Q+22 

0.1008Q+22 
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Table  F-ll. 

Calculated  Solution  of  P2  by  Si 
Stable  Ketncd. 

xth  Order  A{  1 . 

321,  0.121)- 

Step 

First  Component 

Second 

Component 

(h=0.005) 

Real  Imaginary 

Real 

Imaginary 

0 

0.2000Q+01  0.0 

-0.9900Q+QQ 

0. 2500Q*01 

200 

0.3679Q+00  0.5007Q-10 

-0.471 5Q-09 

0.751SQ-C9 

400 

0.1353Q+00  -0.20S9Q-18 

0.4503Q-18 

0.38S5Q-18 

600 

0.4979Q-01  -0.6950Q-28 

0.3034Q-27 

-0.25SGQ-27 

800 

0.1832Q-01  —0.21 82Q-31 

-0.14C8Q-36 

-0.2275Q- 36 

1000 

0.6738Q-02  -0.5375Q-32 

-0 . 1 651 Q— 45 

0.7022Q-46 

Table  F-12. 

Actual  Global  Error  of  Calculated  Solution  of  P2  by  Sixth 
Order  A( 1.321,  C.l2i)-Stable  Method. 

First  Component 

Second  Component 

Step 

Real  Imaginary 

Real 

Imaginary 

.  0 

0.0  0.0 

0.0 

0.0 

200 

-0.3246Q-QS  -0.6007Q-10 

0.471 5Q-09 

-0.7519Q-09 

400 

0.3310Q-14  0.2089Q-18 

-0.4503Q-18 

-0.3885Q-1S 

600 

0.2113Q-14  0. 6950 Q- 28 

-0.3034Q-27 

0.2590Q-27 

800 

0.1039Q-14  0. 2182 Q- 31 

0.1408Q-36 

0.2275Q-36 

1000 

0.4786Q-15  0.5375Q-32 

0.1651  Q-45 

-0.7022Q-46 

1  I 


Table  F-13.  Calculated  Solution  of  P2  by  Sixth  Order  BDF. 


Step 

First 

Component 

Second  Component 

(h=0.005) 

Real 

Imaginary 

Real 

Imaginary 

0 

0.2CGCQ+C1 

0.0 

-0.990QQ+00 

0.2500Q+01 

200 

-0.16S7Q+I4 

-G.431SQ+14 

0.1245Q+15 

0.1 083 Q+l 3 

400 

0.3979Q+2S 

C.5945Q+29 

-0.1 £20  Q+30 

0. 4061 Q+23 

600 

-0.8016Q+44 

-0.7611 Q+44 

0. 2696Q+45 

-0.125GQ+45 

800 

0.1467Q+C0 

0.87530+59 

-0.364GQ+60 

0.2800Q+60 

1000 

-G.25Q8Q+75 

-0.821 6Q+74 

0.4537Q+75 

-0. 5457Q+75 

Table  F-14. 

Actual  Global  Error  of  Calculat 
Order  BDF. 

ed  Solution  of  P2  by  Sixth 

Step 

First  Component 

Second  Component 

Real 

Imaginary 

Peal 

Imaginary 

0 

0.0 

0.0 

0.0 

0.0 

200 

0.1667Q+14 

0.431SQ+14 

-0.1245Q+15 

-0.10890+1 3 

400 

• 

-0.3979Q+29 

-0. 5945Q+29 

*  0.1880Q+30 

-G. 40610+29 

600 

0.8015Q+44 

0.7611Q+44 

-0.2696Q+45 

0.1250Q+45 

800 

-0.1467Q+60 

-0. S753Q+59 

0.3640QJ60 

-O.2S0OQ+6G 

1000 

0.2508Q+75 

0.8216Q+74 

-0.4537Q+75 

0.5457Q+75 

